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

Last change on this file since 1449 was 1446, checked in by TakTsutsumi, 16 years ago

Merged recent updates (since 2007) from nrao-asap
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 90.4 KB
Line 
1//
2// C++ Implementation: STMath
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13#include <casa/iomanip.h>
14#include <casa/Exceptions/Error.h>
15#include <casa/Containers/Block.h>
16#include <casa/BasicSL/String.h>
17#include <casa/Arrays/MaskArrLogi.h>
18#include <casa/Arrays/MaskArrMath.h>
19#include <casa/Arrays/ArrayLogical.h>
20#include <casa/Arrays/ArrayMath.h>
21#include <casa/Arrays/Slice.h>
22#include <casa/Arrays/Slicer.h>
23#include <casa/Containers/RecordField.h>
24#include <tables/Tables/TableRow.h>
25#include <tables/Tables/TableVector.h>
26#include <tables/Tables/TabVecMath.h>
27#include <tables/Tables/ExprNode.h>
28#include <tables/Tables/TableRecord.h>
29#include <tables/Tables/TableParse.h>
30#include <tables/Tables/ReadAsciiTable.h>
31#include <tables/Tables/TableIter.h>
32#include <tables/Tables/TableCopy.h>
33#include <scimath/Mathematics/FFTServer.h>
34
35#include <lattices/Lattices/LatticeUtilities.h>
36
37#include <coordinates/Coordinates/SpectralCoordinate.h>
38#include <coordinates/Coordinates/CoordinateSystem.h>
39#include <coordinates/Coordinates/CoordinateUtil.h>
40#include <coordinates/Coordinates/FrequencyAligner.h>
41
42#include <scimath/Mathematics/VectorKernel.h>
43#include <scimath/Mathematics/Convolver.h>
44#include <scimath/Functionals/Polynomial.h>
45
46#include "MathUtils.h"
47#include "RowAccumulator.h"
48#include "STAttr.h"
49#include "STMath.h"
50#include "STSelector.h"
51
52using namespace casa;
53
54using namespace asap;
55
56STMath::STMath(bool insitu) :
57 insitu_(insitu)
58{
59}
60
61
62STMath::~STMath()
63{
64}
65
66CountedPtr<Scantable>
67STMath::average( const std::vector<CountedPtr<Scantable> >& in,
68 const std::vector<bool>& mask,
69 const std::string& weight,
70 const std::string& avmode)
71{
72 if ( avmode == "SCAN" && in.size() != 1 )
73 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
74 "Use merge first."));
75 WeightType wtype = stringToWeight(weight);
76
77 // output
78 // clone as this is non insitu
79 bool insitu = insitu_;
80 setInsitu(false);
81 CountedPtr< Scantable > out = getScantable(in[0], true);
82 setInsitu(insitu);
83 std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
84 ++stit;
85 while ( stit != in.end() ) {
86 out->appendToHistoryTable((*stit)->history());
87 ++stit;
88 }
89
90 Table& tout = out->table();
91
92 /// @todo check if all scantables are conformant
93
94 ArrayColumn<Float> specColOut(tout,"SPECTRA");
95 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
96 ArrayColumn<Float> tsysColOut(tout,"TSYS");
97 ScalarColumn<Double> mjdColOut(tout,"TIME");
98 ScalarColumn<Double> intColOut(tout,"INTERVAL");
99 ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
100 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
101
102 // set up the output table rows. These are based on the structure of the
103 // FIRST scantable in the vector
104 const Table& baset = in[0]->table();
105
106 Block<String> cols(3);
107 cols[0] = String("BEAMNO");
108 cols[1] = String("IFNO");
109 cols[2] = String("POLNO");
110 if ( avmode == "SOURCE" ) {
111 cols.resize(4);
112 cols[3] = String("SRCNAME");
113 }
114 if ( avmode == "SCAN" && in.size() == 1) {
115 cols.resize(4);
116 cols[3] = String("SCANNO");
117 }
118 uInt outrowCount = 0;
119 TableIterator iter(baset, cols);
120 while (!iter.pastEnd()) {
121 Table subt = iter.table();
122 // copy the first row of this selection into the new table
123 tout.addRow();
124 TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
125 // re-index to 0
126 if ( avmode != "SCAN" && avmode != "SOURCE" ) {
127 scanColOut.put(outrowCount, uInt(0));
128 }
129 ++outrowCount;
130 ++iter;
131 }
132 RowAccumulator acc(wtype);
133 Vector<Bool> cmask(mask);
134 acc.setUserMask(cmask);
135 ROTableRow row(tout);
136 ROArrayColumn<Float> specCol, tsysCol;
137 ROArrayColumn<uChar> flagCol;
138 ROScalarColumn<Double> mjdCol, intCol;
139 ROScalarColumn<Int> scanIDCol;
140
141 Vector<uInt> rowstodelete;
142
143 for (uInt i=0; i < tout.nrow(); ++i) {
144 for ( int j=0; j < int(in.size()); ++j ) {
145 const Table& tin = in[j]->table();
146 const TableRecord& rec = row.get(i);
147 ROScalarColumn<Double> tmp(tin, "TIME");
148 Double td;tmp.get(0,td);
149 Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
150 && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
151 && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
152 Table subt;
153 if ( avmode == "SOURCE") {
154 subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
155 } else if (avmode == "SCAN") {
156 subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
157 } else {
158 subt = basesubt;
159 }
160 specCol.attach(subt,"SPECTRA");
161 flagCol.attach(subt,"FLAGTRA");
162 tsysCol.attach(subt,"TSYS");
163 intCol.attach(subt,"INTERVAL");
164 mjdCol.attach(subt,"TIME");
165 Vector<Float> spec,tsys;
166 Vector<uChar> flag;
167 Double inter,time;
168 for (uInt k = 0; k < subt.nrow(); ++k ) {
169 flagCol.get(k, flag);
170 Vector<Bool> bflag(flag.shape());
171 convertArray(bflag, flag);
172 /*
173 if ( allEQ(bflag, True) ) {
174 continue;//don't accumulate
175 }
176 */
177 specCol.get(k, spec);
178 tsysCol.get(k, tsys);
179 intCol.get(k, inter);
180 mjdCol.get(k, time);
181 // spectrum has to be added last to enable weighting by the other values
182 acc.add(spec, !bflag, tsys, inter, time);
183 }
184 }
185 const Vector<Bool>& msk = acc.getMask();
186 if ( allEQ(msk, False) ) {
187 uint n = rowstodelete.nelements();
188 rowstodelete.resize(n+1, True);
189 rowstodelete[n] = i;
190 continue;
191 }
192 //write out
193 if (acc.state()) {
194 Vector<uChar> flg(msk.shape());
195 convertArray(flg, !msk);
196 flagColOut.put(i, flg);
197 specColOut.put(i, acc.getSpectrum());
198 tsysColOut.put(i, acc.getTsys());
199 intColOut.put(i, acc.getInterval());
200 mjdColOut.put(i, acc.getTime());
201 // we should only have one cycle now -> reset it to be 0
202 // frequency switched data has different CYCLENO for different IFNO
203 // which requires resetting this value
204 cycColOut.put(i, uInt(0));
205 } else {
206 ostringstream oss;
207 oss << "For output row="<<i<<", all input rows of data are flagged. no averaging" << endl;
208 pushLog(String(oss));
209 }
210 acc.reset();
211 }
212 if (rowstodelete.nelements() > 0) {
213 cout << rowstodelete << endl;
214 tout.removeRow(rowstodelete);
215 if (tout.nrow() == 0) {
216 throw(AipsError("Can't average fully flagged data."));
217 }
218 }
219 return out;
220}
221
222CountedPtr< Scantable >
223 STMath::averageChannel( const CountedPtr < Scantable > & in,
224 const std::string & mode,
225 const std::string& avmode )
226{
227 // clone as this is non insitu
228 bool insitu = insitu_;
229 setInsitu(false);
230 CountedPtr< Scantable > out = getScantable(in, true);
231 setInsitu(insitu);
232 Table& tout = out->table();
233 ArrayColumn<Float> specColOut(tout,"SPECTRA");
234 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
235 ArrayColumn<Float> tsysColOut(tout,"TSYS");
236 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
237 ScalarColumn<Double> intColOut(tout, "INTERVAL");
238 Table tmp = in->table().sort("BEAMNO");
239 Block<String> cols(3);
240 cols[0] = String("BEAMNO");
241 cols[1] = String("IFNO");
242 cols[2] = String("POLNO");
243 if ( avmode == "SCAN") {
244 cols.resize(4);
245 cols[3] = String("SCANNO");
246 }
247 uInt outrowCount = 0;
248 uChar userflag = 1 << 7;
249 TableIterator iter(tmp, cols);
250 while (!iter.pastEnd()) {
251 Table subt = iter.table();
252 ROArrayColumn<Float> specCol, tsysCol;
253 ROArrayColumn<uChar> flagCol;
254 ROScalarColumn<Double> intCol(subt, "INTERVAL");
255 specCol.attach(subt,"SPECTRA");
256 flagCol.attach(subt,"FLAGTRA");
257 tsysCol.attach(subt,"TSYS");
258 tout.addRow();
259 TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
260 if ( avmode != "SCAN") {
261 scanColOut.put(outrowCount, uInt(0));
262 }
263 Vector<Float> tmp;
264 specCol.get(0, tmp);
265 uInt nchan = tmp.nelements();
266 // have to do channel by channel here as MaskedArrMath
267 // doesn't have partialMedians
268 Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
269 Vector<Float> outspec(nchan);
270 Vector<uChar> outflag(nchan,0);
271 Vector<Float> outtsys(1);/// @fixme when tsys is channel based
272 for (uInt i=0; i<nchan; ++i) {
273 Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
274 MaskedArray<Float> ma = maskedArray(specs,flags);
275 outspec[i] = median(ma);
276 if ( allEQ(ma.getMask(), False) )
277 outflag[i] = userflag;// flag data
278 }
279 outtsys[0] = median(tsysCol.getColumn());
280 specColOut.put(outrowCount, outspec);
281 flagColOut.put(outrowCount, outflag);
282 tsysColOut.put(outrowCount, outtsys);
283 Double intsum = sum(intCol.getColumn());
284 intColOut.put(outrowCount, intsum);
285 ++outrowCount;
286 ++iter;
287 }
288 return out;
289}
290
291CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
292 bool droprows)
293{
294 if (insitu_) return in;
295 else {
296 // clone
297 Scantable* tabp = new Scantable(*in, Bool(droprows));
298 return CountedPtr<Scantable>(tabp);
299 }
300}
301
302CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
303 float val,
304 const std::string& mode,
305 bool tsys )
306{
307 CountedPtr< Scantable > out = getScantable(in, false);
308 Table& tab = out->table();
309 ArrayColumn<Float> specCol(tab,"SPECTRA");
310 ArrayColumn<Float> tsysCol(tab,"TSYS");
311 for (uInt i=0; i<tab.nrow(); ++i) {
312 Vector<Float> spec;
313 Vector<Float> ts;
314 specCol.get(i, spec);
315 tsysCol.get(i, ts);
316 if (mode == "MUL" || mode == "DIV") {
317 if (mode == "DIV") val = 1.0/val;
318 spec *= val;
319 specCol.put(i, spec);
320 if ( tsys ) {
321 ts *= val;
322 tsysCol.put(i, ts);
323 }
324 } else if ( mode == "ADD" || mode == "SUB") {
325 if (mode == "SUB") val *= -1.0;
326 spec += val;
327 specCol.put(i, spec);
328 if ( tsys ) {
329 ts += val;
330 tsysCol.put(i, ts);
331 }
332 }
333 }
334 return out;
335}
336
337CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
338 const CountedPtr<Scantable>& right,
339 const std::string& mode)
340{
341 bool insitu = insitu_;
342 if ( ! left->conformant(*right) ) {
343 throw(AipsError("'left' and 'right' scantables are not conformant."));
344 }
345 setInsitu(false);
346 CountedPtr< Scantable > out = getScantable(left, false);
347 setInsitu(insitu);
348 Table& tout = out->table();
349 Block<String> coln(5);
350 coln[0] = "SCANNO"; coln[1] = "CYCLENO"; coln[2] = "BEAMNO";
351 coln[3] = "IFNO"; coln[4] = "POLNO";
352 Table tmpl = tout.sort(coln);
353 Table tmpr = right->table().sort(coln);
354 ArrayColumn<Float> lspecCol(tmpl,"SPECTRA");
355 ROArrayColumn<Float> rspecCol(tmpr,"SPECTRA");
356 ArrayColumn<uChar> lflagCol(tmpl,"FLAGTRA");
357 ROArrayColumn<uChar> rflagCol(tmpr,"FLAGTRA");
358
359 for (uInt i=0; i<tout.nrow(); ++i) {
360 Vector<Float> lspecvec, rspecvec;
361 Vector<uChar> lflagvec, rflagvec;
362 lspecvec = lspecCol(i); rspecvec = rspecCol(i);
363 lflagvec = lflagCol(i); rflagvec = rflagCol(i);
364 MaskedArray<Float> mleft = maskedArray(lspecvec, lflagvec);
365 MaskedArray<Float> mright = maskedArray(rspecvec, rflagvec);
366 if (mode == "ADD") {
367 mleft += mright;
368 } else if ( mode == "SUB") {
369 mleft -= mright;
370 } else if ( mode == "MUL") {
371 mleft *= mright;
372 } else if ( mode == "DIV") {
373 mleft /= mright;
374 } else {
375 throw(AipsError("Illegal binary operator"));
376 }
377 lspecCol.put(i, mleft.getArray());
378 }
379 return out;
380}
381
382
383
384MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
385 const Vector<uChar>& f)
386{
387 Vector<Bool> mask;
388 mask.resize(f.shape());
389 convertArray(mask, f);
390 return MaskedArray<Float>(s,!mask);
391}
392
393Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
394{
395 const Vector<Bool>& m = ma.getMask();
396 Vector<uChar> flags(m.shape());
397 convertArray(flags, !m);
398 return flags;
399}
400
401CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
402 const std::string & mode,
403 bool preserve )
404{
405 /// @todo make other modes available
406 /// modes should be "nearest", "pair"
407 // make this operation non insitu
408 const Table& tin = in->table();
409 Table ons = tin(tin.col("SRCTYPE") == Int(0));
410 Table offs = tin(tin.col("SRCTYPE") == Int(1));
411 if ( offs.nrow() == 0 )
412 throw(AipsError("No 'off' scans present."));
413 // put all "on" scans into output table
414
415 bool insitu = insitu_;
416 setInsitu(false);
417 CountedPtr< Scantable > out = getScantable(in, true);
418 setInsitu(insitu);
419 Table& tout = out->table();
420
421 TableCopy::copyRows(tout, ons);
422 TableRow row(tout);
423 ROScalarColumn<Double> offtimeCol(offs, "TIME");
424 ArrayColumn<Float> outspecCol(tout, "SPECTRA");
425 ROArrayColumn<Float> outtsysCol(tout, "TSYS");
426 ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
427 for (uInt i=0; i < tout.nrow(); ++i) {
428 const TableRecord& rec = row.get(i);
429 Double ontime = rec.asDouble("TIME");
430 Table presel = offs(offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
431 && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
432 && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
433 ROScalarColumn<Double> offtimeCol(presel, "TIME");
434
435 Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
436 // Timestamp may vary within a cycle ???!!!
437 // increase this by 0.01 sec in case of rounding errors...
438 // There might be a better way to do this.
439 // fix to this fix. TIME is MJD, so 1.0d not 1.0s
440 mindeltat += 0.01/24./60./60.;
441 Table sel = presel( abs(presel.col("TIME")-ontime) <= mindeltat);
442
443 if ( sel.nrow() < 1 ) {
444 throw(AipsError("No closest in time found... This could be a rounding "
445 "issue. Try quotient instead."));
446 }
447 TableRow offrow(sel);
448 const TableRecord& offrec = offrow.get(0);//should only be one row
449 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
450 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
451 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
452 /// @fixme this assumes tsys is a scalar not vector
453 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
454 Vector<Float> specon, tsyson;
455 outtsysCol.get(i, tsyson);
456 outspecCol.get(i, specon);
457 Vector<uChar> flagon;
458 outflagCol.get(i, flagon);
459 MaskedArray<Float> mon = maskedArray(specon, flagon);
460 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
461 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
462 if (preserve) {
463 quot -= tsysoffscalar;
464 } else {
465 quot -= tsyson[0];
466 }
467 outspecCol.put(i, quot.getArray());
468 outflagCol.put(i, flagsFromMA(quot));
469 }
470 // renumber scanno
471 TableIterator it(tout, "SCANNO");
472 uInt i = 0;
473 while ( !it.pastEnd() ) {
474 Table t = it.table();
475 TableVector<uInt> vec(t, "SCANNO");
476 vec = i;
477 ++i;
478 ++it;
479 }
480 return out;
481}
482
483
484CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
485 const CountedPtr< Scantable > & off,
486 bool preserve )
487{
488 bool insitu = insitu_;
489 if ( ! on->conformant(*off) ) {
490 throw(AipsError("'on' and 'off' scantables are not conformant."));
491 }
492 setInsitu(false);
493 CountedPtr< Scantable > out = getScantable(on, false);
494 setInsitu(insitu);
495 Table& tout = out->table();
496 const Table& toff = off->table();
497 TableIterator sit(tout, "SCANNO");
498 TableIterator s2it(toff, "SCANNO");
499 while ( !sit.pastEnd() ) {
500 Table ton = sit.table();
501 TableRow row(ton);
502 Table t = s2it.table();
503 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
504 ROArrayColumn<Float> outtsysCol(ton, "TSYS");
505 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
506 for (uInt i=0; i < ton.nrow(); ++i) {
507 const TableRecord& rec = row.get(i);
508 Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
509 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
510 && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
511 if ( offsel.nrow() == 0 )
512 throw AipsError("STMath::quotient: no matching off");
513 TableRow offrow(offsel);
514 const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
515 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
516 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
517 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
518 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
519 Vector<Float> specon, tsyson;
520 outtsysCol.get(i, tsyson);
521 outspecCol.get(i, specon);
522 Vector<uChar> flagon;
523 outflagCol.get(i, flagon);
524 MaskedArray<Float> mon = maskedArray(specon, flagon);
525 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
526 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
527 if (preserve) {
528 quot -= tsysoffscalar;
529 } else {
530 quot -= tsyson[0];
531 }
532 outspecCol.put(i, quot.getArray());
533 outflagCol.put(i, flagsFromMA(quot));
534 }
535 ++sit;
536 ++s2it;
537 // take the first off for each on scan which doesn't have a
538 // matching off scan
539 // non <= noff: matching pairs, non > noff matching pairs then first off
540 if ( s2it.pastEnd() ) s2it.reset();
541 }
542 return out;
543}
544
545// dototalpower (migration of GBTIDL procedure dototalpower.pro)
546// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
547// do it for each cycles in a specific scan.
548CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
549 const CountedPtr< Scantable >& caloff, Float tcal )
550{
551if ( ! calon->conformant(*caloff) ) {
552 throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
553 }
554 setInsitu(false);
555 CountedPtr< Scantable > out = getScantable(caloff, false);
556 Table& tout = out->table();
557 const Table& tcon = calon->table();
558 Vector<Float> tcalout;
559 Vector<Float> tcalout2; //debug
560
561 if ( tout.nrow() != tcon.nrow() ) {
562 throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
563 }
564 // iteration by scanno or cycle no.
565 TableIterator sit(tout, "SCANNO");
566 TableIterator s2it(tcon, "SCANNO");
567 while ( !sit.pastEnd() ) {
568 Table toff = sit.table();
569 TableRow row(toff);
570 Table t = s2it.table();
571 ScalarColumn<Double> outintCol(toff, "INTERVAL");
572 ArrayColumn<Float> outspecCol(toff, "SPECTRA");
573 ArrayColumn<Float> outtsysCol(toff, "TSYS");
574 ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
575 ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
576 ROScalarColumn<uInt> outpolCol(toff, "POLNO");
577 ROScalarColumn<Double> onintCol(t, "INTERVAL");
578 ROArrayColumn<Float> onspecCol(t, "SPECTRA");
579 ROArrayColumn<Float> ontsysCol(t, "TSYS");
580 ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
581 //ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
582
583 for (uInt i=0; i < toff.nrow(); ++i) {
584 //skip these checks -> assumes the data order are the same between the cal on off pairs
585 //
586 Vector<Float> specCalon, specCaloff;
587 // to store scalar (mean) tsys
588 Vector<Float> tsysout(1);
589 uInt tcalId, polno;
590 Double offint, onint;
591 outpolCol.get(i, polno);
592 outspecCol.get(i, specCaloff);
593 onspecCol.get(i, specCalon);
594 Vector<uChar> flagCaloff, flagCalon;
595 outflagCol.get(i, flagCaloff);
596 onflagCol.get(i, flagCalon);
597 outtcalIdCol.get(i, tcalId);
598 outintCol.get(i, offint);
599 onintCol.get(i, onint);
600 // caluculate mean Tsys
601 uInt nchan = specCaloff.nelements();
602 // percentage of edge cut off
603 uInt pc = 10;
604 uInt bchan = nchan/pc;
605 uInt echan = nchan-bchan;
606
607 Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
608 Vector<Float> testsubsp = specCaloff(chansl);
609 MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
610 MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
611 MaskedArray<Float> spdiff = spon-spoff;
612 uInt noff = spoff.nelementsValid();
613 //uInt non = spon.nelementsValid();
614 uInt ndiff = spdiff.nelementsValid();
615 Float meantsys;
616
617/**
618 Double subspec, subdiff;
619 uInt usednchan;
620 subspec = 0;
621 subdiff = 0;
622 usednchan = 0;
623 for(uInt k=(bchan-1); k<echan; k++) {
624 subspec += specCaloff[k];
625 subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
626 ++usednchan;
627 }
628**/
629 // get tcal if input tcal <= 0
630 String tcalt;
631 Float tcalUsed;
632 tcalUsed = tcal;
633 if ( tcal <= 0.0 ) {
634 caloff->tcal().getEntry(tcalt, tcalout, tcalId);
635 if (polno<=3) {
636 tcalUsed = tcalout[polno];
637 }
638 else {
639 tcalUsed = tcalout[0];
640 }
641 }
642
643 Float meanoff;
644 Float meandiff;
645 if (noff && ndiff) {
646 //Debug
647 //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
648 meanoff = sum(spoff)/noff;
649 meandiff = sum(spdiff)/ndiff;
650 meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
651 }
652 else {
653 meantsys=1;
654 }
655
656 tsysout[0] = Float(meantsys);
657 MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
658 MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
659 MaskedArray<Float> sig = Float(0.5) * (mcaloff + mcalon);
660 //uInt ncaloff = mcaloff.nelementsValid();
661 //uInt ncalon = mcalon.nelementsValid();
662
663 outintCol.put(i, offint+onint);
664 outspecCol.put(i, sig.getArray());
665 outflagCol.put(i, flagsFromMA(sig));
666 outtsysCol.put(i, tsysout);
667 }
668 ++sit;
669 ++s2it;
670 }
671 return out;
672}
673
674//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
675// observatiions.
676// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
677// no smoothing).
678// output: resultant scantable [= (sig-ref/ref)*tsys]
679CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
680 const CountedPtr < Scantable >& ref,
681 int smoothref,
682 casa::Float tsysv,
683 casa::Float tau )
684{
685if ( ! ref->conformant(*sig) ) {
686 throw(AipsError("'sig' and 'ref' scantables are not conformant."));
687 }
688 setInsitu(false);
689 CountedPtr< Scantable > out = getScantable(sig, false);
690 CountedPtr< Scantable > smref;
691 if ( smoothref > 1 ) {
692 float fsmoothref = static_cast<float>(smoothref);
693 std::string inkernel = "boxcar";
694 smref = smooth(ref, inkernel, fsmoothref );
695 ostringstream oss;
696 oss<<"Applied smoothing of "<<fsmoothref<<" on the reference."<<endl;
697 pushLog(String(oss));
698 }
699 else {
700 smref = ref;
701 }
702 Table& tout = out->table();
703 const Table& tref = smref->table();
704 if ( tout.nrow() != tref.nrow() ) {
705 throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
706 }
707 // iteration by scanno? or cycle no.
708 TableIterator sit(tout, "SCANNO");
709 TableIterator s2it(tref, "SCANNO");
710 while ( !sit.pastEnd() ) {
711 Table ton = sit.table();
712 Table t = s2it.table();
713 ScalarColumn<Double> outintCol(ton, "INTERVAL");
714 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
715 ArrayColumn<Float> outtsysCol(ton, "TSYS");
716 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
717 ArrayColumn<Float> refspecCol(t, "SPECTRA");
718 ROScalarColumn<Double> refintCol(t, "INTERVAL");
719 ROArrayColumn<Float> reftsysCol(t, "TSYS");
720 ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
721 ROScalarColumn<Float> refelevCol(t, "ELEVATION");
722 for (uInt i=0; i < ton.nrow(); ++i) {
723
724 Double onint, refint;
725 Vector<Float> specon, specref;
726 // to store scalar (mean) tsys
727 Vector<Float> tsysref;
728 outintCol.get(i, onint);
729 refintCol.get(i, refint);
730 outspecCol.get(i, specon);
731 refspecCol.get(i, specref);
732 Vector<uChar> flagref, flagon;
733 outflagCol.get(i, flagon);
734 refflagCol.get(i, flagref);
735 reftsysCol.get(i, tsysref);
736
737 Float tsysrefscalar;
738 if ( tsysv > 0.0 ) {
739 ostringstream oss;
740 Float elev;
741 refelevCol.get(i, elev);
742 oss << "user specified Tsys = " << tsysv;
743 // do recalc elevation if EL = 0
744 if ( elev == 0 ) {
745 throw(AipsError("EL=0, elevation data is missing."));
746 } else {
747 if ( tau <= 0.0 ) {
748 throw(AipsError("Valid tau is not supplied."));
749 } else {
750 tsysrefscalar = tsysv * exp(tau/elev);
751 }
752 }
753 oss << ", corrected (for El) tsys= "<<tsysrefscalar;
754 pushLog(String(oss));
755 }
756 else {
757 tsysrefscalar = tsysref[0];
758 }
759 //get quotient spectrum
760 MaskedArray<Float> mref = maskedArray(specref, flagref);
761 MaskedArray<Float> mon = maskedArray(specon, flagon);
762 MaskedArray<Float> specres = tsysrefscalar*((mon - mref)/mref);
763 Double resint = onint*refint*smoothref/(onint+refint*smoothref);
764
765 //Debug
766 //cerr<<"Tsys used="<<tsysrefscalar<<endl;
767 // fill the result, replay signal tsys by reference tsys
768 outintCol.put(i, resint);
769 outspecCol.put(i, specres.getArray());
770 outflagCol.put(i, flagsFromMA(specres));
771 outtsysCol.put(i, tsysref);
772 }
773 ++sit;
774 ++s2it;
775 }
776 return out;
777}
778
779CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
780 const std::vector<int>& scans,
781 int smoothref,
782 casa::Float tsysv,
783 casa::Float tau,
784 casa::Float tcal )
785
786{
787 setInsitu(false);
788 STSelector sel;
789 std::vector<int> scan1, scan2, beams;
790 std::vector< vector<int> > scanpair;
791 std::vector<string> calstate;
792 String msg;
793
794 CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
795 CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
796
797 std::vector< CountedPtr< Scantable > > sctables;
798 sctables.push_back(s1b1on);
799 sctables.push_back(s1b1off);
800 sctables.push_back(s1b2on);
801 sctables.push_back(s1b2off);
802 sctables.push_back(s2b1on);
803 sctables.push_back(s2b1off);
804 sctables.push_back(s2b2on);
805 sctables.push_back(s2b2off);
806
807 //check scanlist
808 int n=s->checkScanInfo(scans);
809 if (n==1) {
810 throw(AipsError("Incorrect scan pairs. "));
811 }
812
813 // Assume scans contain only a pair of consecutive scan numbers.
814 // It is assumed that first beam, b1, is on target.
815 // There is no check if the first beam is on or not.
816 if ( scans.size()==1 ) {
817 scan1.push_back(scans[0]);
818 scan2.push_back(scans[0]+1);
819 } else if ( scans.size()==2 ) {
820 scan1.push_back(scans[0]);
821 scan2.push_back(scans[1]);
822 } else {
823 if ( scans.size()%2 == 0 ) {
824 for (uInt i=0; i<scans.size(); i++) {
825 if (i%2 == 0) {
826 scan1.push_back(scans[i]);
827 }
828 else {
829 scan2.push_back(scans[i]);
830 }
831 }
832 } else {
833 throw(AipsError("Odd numbers of scans, cannot form pairs."));
834 }
835 }
836 scanpair.push_back(scan1);
837 scanpair.push_back(scan2);
838 calstate.push_back("*calon");
839 calstate.push_back("*[^calon]");
840 CountedPtr< Scantable > ws = getScantable(s, false);
841 uInt l=0;
842 while ( l < sctables.size() ) {
843 for (uInt i=0; i < 2; i++) {
844 for (uInt j=0; j < 2; j++) {
845 for (uInt k=0; k < 2; k++) {
846 sel.reset();
847 sel.setScans(scanpair[i]);
848 sel.setName(calstate[k]);
849 beams.clear();
850 beams.push_back(j);
851 sel.setBeams(beams);
852 ws->setSelection(sel);
853 sctables[l]= getScantable(ws, false);
854 l++;
855 }
856 }
857 }
858 }
859
860 // replace here by splitData or getData functionality
861 CountedPtr< Scantable > sig1;
862 CountedPtr< Scantable > ref1;
863 CountedPtr< Scantable > sig2;
864 CountedPtr< Scantable > ref2;
865 CountedPtr< Scantable > calb1;
866 CountedPtr< Scantable > calb2;
867
868 msg=String("Processing dototalpower for subset of the data");
869 ostringstream oss1;
870 oss1 << msg << endl;
871 pushLog(String(oss1));
872 // Debug for IRC CS data
873 //float tcal1=7.0;
874 //float tcal2=4.0;
875 sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
876 ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
877 ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
878 sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
879
880 // correction of user-specified tsys for elevation here
881
882 // dosigref calibration
883 msg=String("Processing dosigref for subset of the data");
884 ostringstream oss2;
885 oss2 << msg << endl;
886 pushLog(String(oss2));
887 calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
888 calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
889
890 // iteration by scanno or cycle no.
891 Table& tcalb1 = calb1->table();
892 Table& tcalb2 = calb2->table();
893 TableIterator sit(tcalb1, "SCANNO");
894 TableIterator s2it(tcalb2, "SCANNO");
895 while ( !sit.pastEnd() ) {
896 Table t1 = sit.table();
897 Table t2= s2it.table();
898 ArrayColumn<Float> outspecCol(t1, "SPECTRA");
899 ArrayColumn<Float> outtsysCol(t1, "TSYS");
900 ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
901 ScalarColumn<Double> outintCol(t1, "INTERVAL");
902 ArrayColumn<Float> t2specCol(t2, "SPECTRA");
903 ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
904 ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
905 ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
906 for (uInt i=0; i < t1.nrow(); ++i) {
907 Vector<Float> spec1, spec2;
908 // to store scalar (mean) tsys
909 Vector<Float> tsys1, tsys2;
910 Vector<uChar> flag1, flag2;
911 Double tint1, tint2;
912 outspecCol.get(i, spec1);
913 t2specCol.get(i, spec2);
914 outflagCol.get(i, flag1);
915 t2flagCol.get(i, flag2);
916 outtsysCol.get(i, tsys1);
917 t2tsysCol.get(i, tsys2);
918 outintCol.get(i, tint1);
919 t2intCol.get(i, tint2);
920 // average
921 // assume scalar tsys for weights
922 Float wt1, wt2, tsyssq1, tsyssq2;
923 tsyssq1 = tsys1[0]*tsys1[0];
924 tsyssq2 = tsys2[0]*tsys2[0];
925 wt1 = Float(tint1)/tsyssq1;
926 wt2 = Float(tint2)/tsyssq2;
927 Float invsumwt=1/(wt1+wt2);
928 MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
929 MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
930 MaskedArray<Float> avspec = invsumwt * (wt1*mspec1 + wt2*mspec2);
931 //Array<Float> avtsys = Float(0.5) * (tsys1 + tsys2);
932 // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
933 tsys1[0] = sqrt(tsyssq1 + tsyssq2);
934 Array<Float> avtsys = tsys1;
935
936 outspecCol.put(i, avspec.getArray());
937 outflagCol.put(i, flagsFromMA(avspec));
938 outtsysCol.put(i, avtsys);
939 }
940 ++sit;
941 ++s2it;
942 }
943 return calb1;
944}
945
946//GBTIDL version of frequency switched data calibration
947CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
948 const std::vector<int>& scans,
949 int smoothref,
950 casa::Float tsysv,
951 casa::Float tau,
952 casa::Float tcal )
953{
954
955
956 STSelector sel;
957 CountedPtr< Scantable > ws = getScantable(s, false);
958 CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
959 CountedPtr< Scantable > calsig, calref, out, out1, out2;
960 Bool nofold=False;
961
962 //split the data
963 sel.setName("*_fs");
964 ws->setSelection(sel);
965 sig = getScantable(ws,false);
966 sel.reset();
967 sel.setName("*_fs_calon");
968 ws->setSelection(sel);
969 sigwcal = getScantable(ws,false);
970 sel.reset();
971 sel.setName("*_fsr");
972 ws->setSelection(sel);
973 ref = getScantable(ws,false);
974 sel.reset();
975 sel.setName("*_fsr_calon");
976 ws->setSelection(sel);
977 refwcal = getScantable(ws,false);
978
979 calsig = dototalpower(sigwcal, sig, tcal=tcal);
980 calref = dototalpower(refwcal, ref, tcal=tcal);
981
982 out1=dosigref(calsig,calref,smoothref,tsysv,tau);
983 out2=dosigref(calref,calsig,smoothref,tsysv,tau);
984
985 Table& tabout1=out1->table();
986 Table& tabout2=out2->table();
987 ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
988 ROScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
989 ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
990 Vector<Float> spec; specCol.get(0, spec);
991 uInt nchan = spec.nelements();
992 uInt freqid1; freqidCol1.get(0,freqid1);
993 uInt freqid2; freqidCol2.get(0,freqid2);
994 Double rp1, rp2, rv1, rv2, inc1, inc2;
995 out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
996 out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
997 if (rp1==rp2) {
998 Double foffset = rv1 - rv2;
999 uInt choffset = static_cast<uInt>(foffset/abs(inc2));
1000 if (choffset >= nchan) {
1001 cerr<<"out-band frequency switching, no folding"<<endl;
1002 nofold = True;
1003 }
1004 }
1005
1006 if (nofold) {
1007 std::vector< CountedPtr< Scantable > > tabs;
1008 tabs.push_back(out1);
1009 tabs.push_back(out2);
1010 out = merge(tabs);
1011 }
1012 else { //folding is not implemented yet
1013 out = out1;
1014 }
1015
1016 return out;
1017}
1018
1019CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
1020{
1021 // make copy or reference
1022 CountedPtr< Scantable > out = getScantable(in, false);
1023 Table& tout = out->table();
1024 Block<String> cols(4);
1025 cols[0] = String("SCANNO");
1026 cols[1] = String("CYCLENO");
1027 cols[2] = String("BEAMNO");
1028 cols[3] = String("POLNO");
1029 TableIterator iter(tout, cols);
1030 while (!iter.pastEnd()) {
1031 Table subt = iter.table();
1032 // this should leave us with two rows for the two IFs....if not ignore
1033 if (subt.nrow() != 2 ) {
1034 continue;
1035 }
1036 ArrayColumn<Float> specCol(subt, "SPECTRA");
1037 ArrayColumn<Float> tsysCol(subt, "TSYS");
1038 ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
1039 Vector<Float> onspec,offspec, ontsys, offtsys;
1040 Vector<uChar> onflag, offflag;
1041 tsysCol.get(0, ontsys); tsysCol.get(1, offtsys);
1042 specCol.get(0, onspec); specCol.get(1, offspec);
1043 flagCol.get(0, onflag); flagCol.get(1, offflag);
1044 MaskedArray<Float> on = maskedArray(onspec, onflag);
1045 MaskedArray<Float> off = maskedArray(offspec, offflag);
1046 MaskedArray<Float> oncopy = on.copy();
1047
1048 on /= off; on -= 1.0f;
1049 on *= ontsys[0];
1050 off /= oncopy; off -= 1.0f;
1051 off *= offtsys[0];
1052 specCol.put(0, on.getArray());
1053 const Vector<Bool>& m0 = on.getMask();
1054 Vector<uChar> flags0(m0.shape());
1055 convertArray(flags0, !m0);
1056 flagCol.put(0, flags0);
1057
1058 specCol.put(1, off.getArray());
1059 const Vector<Bool>& m1 = off.getMask();
1060 Vector<uChar> flags1(m1.shape());
1061 convertArray(flags1, !m1);
1062 flagCol.put(1, flags1);
1063 ++iter;
1064 }
1065
1066 return out;
1067}
1068
1069std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
1070 const std::vector< bool > & mask,
1071 const std::string& which )
1072{
1073
1074 Vector<Bool> m(mask);
1075 const Table& tab = in->table();
1076 ROArrayColumn<Float> specCol(tab, "SPECTRA");
1077 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1078 std::vector<float> out;
1079 for (uInt i=0; i < tab.nrow(); ++i ) {
1080 Vector<Float> spec; specCol.get(i, spec);
1081 Vector<uChar> flag; flagCol.get(i, flag);
1082 MaskedArray<Float> ma = maskedArray(spec, flag);
1083 float outstat = 0.0;
1084 if ( spec.nelements() == m.nelements() ) {
1085 outstat = mathutil::statistics(which, ma(m));
1086 } else {
1087 outstat = mathutil::statistics(which, ma);
1088 }
1089 out.push_back(outstat);
1090 }
1091 return out;
1092}
1093
1094CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
1095 int width )
1096{
1097 if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
1098 CountedPtr< Scantable > out = getScantable(in, false);
1099 Table& tout = out->table();
1100 out->frequencies().rescale(width, "BIN");
1101 ArrayColumn<Float> specCol(tout, "SPECTRA");
1102 ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
1103 for (uInt i=0; i < tout.nrow(); ++i ) {
1104 MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
1105 MaskedArray<Float> maout;
1106 LatticeUtilities::bin(maout, main, 0, Int(width));
1107 /// @todo implement channel based tsys binning
1108 specCol.put(i, maout.getArray());
1109 flagCol.put(i, flagsFromMA(maout));
1110 // take only the first binned spectrum's length for the deprecated
1111 // global header item nChan
1112 if (i==0) tout.rwKeywordSet().define(String("nChan"),
1113 Int(maout.getArray().nelements()));
1114 }
1115 return out;
1116}
1117
1118CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
1119 const std::string& method,
1120 float width )
1121//
1122// Should add the possibility of width being specified in km/s. This means
1123// that for each freqID (SpectralCoordinate) we will need to convert to an
1124// average channel width (say at the reference pixel). Then we would need
1125// to be careful to make sure each spectrum (of different freqID)
1126// is the same length.
1127//
1128{
1129 //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1130 Int interpMethod(stringToIMethod(method));
1131
1132 CountedPtr< Scantable > out = getScantable(in, false);
1133 Table& tout = out->table();
1134
1135// Resample SpectralCoordinates (one per freqID)
1136 out->frequencies().rescale(width, "RESAMPLE");
1137 TableIterator iter(tout, "IFNO");
1138 TableRow row(tout);
1139 while ( !iter.pastEnd() ) {
1140 Table tab = iter.table();
1141 ArrayColumn<Float> specCol(tab, "SPECTRA");
1142 //ArrayColumn<Float> tsysCol(tout, "TSYS");
1143 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1144 Vector<Float> spec;
1145 Vector<uChar> flag;
1146 specCol.get(0,spec); // the number of channels should be constant per IF
1147 uInt nChanIn = spec.nelements();
1148 Vector<Float> xIn(nChanIn); indgen(xIn);
1149 Int fac = Int(nChanIn/width);
1150 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
1151 uInt k = 0;
1152 Float x = 0.0;
1153 while (x < Float(nChanIn) ) {
1154 xOut(k) = x;
1155 k++;
1156 x += width;
1157 }
1158 uInt nChanOut = k;
1159 xOut.resize(nChanOut, True);
1160 // process all rows for this IFNO
1161 Vector<Float> specOut;
1162 Vector<Bool> maskOut;
1163 Vector<uChar> flagOut;
1164 for (uInt i=0; i < tab.nrow(); ++i) {
1165 specCol.get(i, spec);
1166 flagCol.get(i, flag);
1167 Vector<Bool> mask(flag.nelements());
1168 convertArray(mask, flag);
1169
1170 IPosition shapeIn(spec.shape());
1171 //sh.nchan = nChanOut;
1172 InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
1173 xIn, spec, mask,
1174 interpMethod, True, True);
1175 /// @todo do the same for channel based Tsys
1176 flagOut.resize(maskOut.nelements());
1177 convertArray(flagOut, maskOut);
1178 specCol.put(i, specOut);
1179 flagCol.put(i, flagOut);
1180 }
1181 ++iter;
1182 }
1183
1184 return out;
1185}
1186
1187STMath::imethod STMath::stringToIMethod(const std::string& in)
1188{
1189 static STMath::imap lookup;
1190
1191 // initialize the lookup table if necessary
1192 if ( lookup.empty() ) {
1193 lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
1194 lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
1195 lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
1196 lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
1197 }
1198
1199 STMath::imap::const_iterator iter = lookup.find(in);
1200
1201 if ( lookup.end() == iter ) {
1202 std::string message = in;
1203 message += " is not a valid interpolation mode";
1204 throw(AipsError(message));
1205 }
1206 return iter->second;
1207}
1208
1209WeightType STMath::stringToWeight(const std::string& in)
1210{
1211 static std::map<std::string, WeightType> lookup;
1212
1213 // initialize the lookup table if necessary
1214 if ( lookup.empty() ) {
1215 lookup["NONE"] = asap::NONE;
1216 lookup["TINT"] = asap::TINT;
1217 lookup["TINTSYS"] = asap::TINTSYS;
1218 lookup["TSYS"] = asap::TSYS;
1219 lookup["VAR"] = asap::VAR;
1220 }
1221
1222 std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
1223
1224 if ( lookup.end() == iter ) {
1225 std::string message = in;
1226 message += " is not a valid weighting mode";
1227 throw(AipsError(message));
1228 }
1229 return iter->second;
1230}
1231
1232CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
1233 const vector< float > & coeff,
1234 const std::string & filename,
1235 const std::string& method)
1236{
1237 // Get elevation data from Scantable and convert to degrees
1238 CountedPtr< Scantable > out = getScantable(in, false);
1239 Table& tab = out->table();
1240 ROScalarColumn<Float> elev(tab, "ELEVATION");
1241 Vector<Float> x = elev.getColumn();
1242 x *= Float(180 / C::pi); // Degrees
1243
1244 Vector<Float> coeffs(coeff);
1245 const uInt nc = coeffs.nelements();
1246 if ( filename.length() > 0 && nc > 0 ) {
1247 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
1248 }
1249
1250 // Correct
1251 if ( nc > 0 || filename.length() == 0 ) {
1252 // Find instrument
1253 Bool throwit = True;
1254 Instrument inst =
1255 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1256 throwit);
1257
1258 // Set polynomial
1259 Polynomial<Float>* ppoly = 0;
1260 Vector<Float> coeff;
1261 String msg;
1262 if ( nc > 0 ) {
1263 ppoly = new Polynomial<Float>(nc);
1264 coeff = coeffs;
1265 msg = String("user");
1266 } else {
1267 STAttr sdAttr;
1268 coeff = sdAttr.gainElevationPoly(inst);
1269 ppoly = new Polynomial<Float>(3);
1270 msg = String("built in");
1271 }
1272
1273 if ( coeff.nelements() > 0 ) {
1274 ppoly->setCoefficients(coeff);
1275 } else {
1276 delete ppoly;
1277 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
1278 }
1279 ostringstream oss;
1280 oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
1281 oss << " " << coeff;
1282 pushLog(String(oss));
1283 const uInt nrow = tab.nrow();
1284 Vector<Float> factor(nrow);
1285 for ( uInt i=0; i < nrow; ++i ) {
1286 factor[i] = 1.0 / (*ppoly)(x[i]);
1287 }
1288 delete ppoly;
1289 scaleByVector(tab, factor, true);
1290
1291 } else {
1292 // Read and correct
1293 pushLog("Making correction from ascii Table");
1294 scaleFromAsciiTable(tab, filename, method, x, true);
1295 }
1296 return out;
1297}
1298
1299void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
1300 const std::string& method,
1301 const Vector<Float>& xout, bool dotsys)
1302{
1303
1304// Read gain-elevation ascii file data into a Table.
1305
1306 String formatString;
1307 Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
1308 scaleFromTable(in, tbl, method, xout, dotsys);
1309}
1310
1311void STMath::scaleFromTable(Table& in,
1312 const Table& table,
1313 const std::string& method,
1314 const Vector<Float>& xout, bool dotsys)
1315{
1316
1317 ROScalarColumn<Float> geElCol(table, "ELEVATION");
1318 ROScalarColumn<Float> geFacCol(table, "FACTOR");
1319 Vector<Float> xin = geElCol.getColumn();
1320 Vector<Float> yin = geFacCol.getColumn();
1321 Vector<Bool> maskin(xin.nelements(),True);
1322
1323 // Interpolate (and extrapolate) with desired method
1324
1325 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1326
1327 Vector<Float> yout;
1328 Vector<Bool> maskout;
1329 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
1330 xin, yin, maskin, interp,
1331 True, True);
1332
1333 scaleByVector(in, Float(1.0)/yout, dotsys);
1334}
1335
1336void STMath::scaleByVector( Table& in,
1337 const Vector< Float >& factor,
1338 bool dotsys )
1339{
1340 uInt nrow = in.nrow();
1341 if ( factor.nelements() != nrow ) {
1342 throw(AipsError("factors.nelements() != table.nelements()"));
1343 }
1344 ArrayColumn<Float> specCol(in, "SPECTRA");
1345 ArrayColumn<uChar> flagCol(in, "FLAGTRA");
1346 ArrayColumn<Float> tsysCol(in, "TSYS");
1347 for (uInt i=0; i < nrow; ++i) {
1348 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1349 ma *= factor[i];
1350 specCol.put(i, ma.getArray());
1351 flagCol.put(i, flagsFromMA(ma));
1352 if ( dotsys ) {
1353 Vector<Float> tsys = tsysCol(i);
1354 tsys *= factor[i];
1355 tsysCol.put(i,tsys);
1356 }
1357 }
1358}
1359
1360CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
1361 float d, float etaap,
1362 float jyperk )
1363{
1364 CountedPtr< Scantable > out = getScantable(in, false);
1365 Table& tab = in->table();
1366 Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
1367 Unit K(String("K"));
1368 Unit JY(String("Jy"));
1369
1370 bool tokelvin = true;
1371 Double cfac = 1.0;
1372
1373 if ( fluxUnit == JY ) {
1374 pushLog("Converting to K");
1375 Quantum<Double> t(1.0,fluxUnit);
1376 Quantum<Double> t2 = t.get(JY);
1377 cfac = (t2 / t).getValue(); // value to Jy
1378
1379 tokelvin = true;
1380 out->setFluxUnit("K");
1381 } else if ( fluxUnit == K ) {
1382 pushLog("Converting to Jy");
1383 Quantum<Double> t(1.0,fluxUnit);
1384 Quantum<Double> t2 = t.get(K);
1385 cfac = (t2 / t).getValue(); // value to K
1386
1387 tokelvin = false;
1388 out->setFluxUnit("Jy");
1389 } else {
1390 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
1391 }
1392 // Make sure input values are converted to either Jy or K first...
1393 Float factor = cfac;
1394
1395 // Select method
1396 if (jyperk > 0.0) {
1397 factor *= jyperk;
1398 if ( tokelvin ) factor = 1.0 / jyperk;
1399 ostringstream oss;
1400 oss << "Jy/K = " << jyperk;
1401 pushLog(String(oss));
1402 Vector<Float> factors(tab.nrow(), factor);
1403 scaleByVector(tab,factors, false);
1404 } else if ( etaap > 0.0) {
1405 if (d < 0) {
1406 Instrument inst =
1407 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1408 True);
1409 STAttr sda;
1410 d = sda.diameter(inst);
1411 }
1412 jyperk = STAttr::findJyPerK(etaap, d);
1413 ostringstream oss;
1414 oss << "Jy/K = " << jyperk;
1415 pushLog(String(oss));
1416 factor *= jyperk;
1417 if ( tokelvin ) {
1418 factor = 1.0 / factor;
1419 }
1420 Vector<Float> factors(tab.nrow(), factor);
1421 scaleByVector(tab, factors, False);
1422 } else {
1423
1424 // OK now we must deal with automatic look up of values.
1425 // We must also deal with the fact that the factors need
1426 // to be computed per IF and may be different and may
1427 // change per integration.
1428
1429 pushLog("Looking up conversion factors");
1430 convertBrightnessUnits(out, tokelvin, cfac);
1431 }
1432
1433 return out;
1434}
1435
1436void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
1437 bool tokelvin, float cfac )
1438{
1439 Table& table = in->table();
1440 Instrument inst =
1441 STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
1442 TableIterator iter(table, "FREQ_ID");
1443 STFrequencies stfreqs = in->frequencies();
1444 STAttr sdAtt;
1445 while (!iter.pastEnd()) {
1446 Table tab = iter.table();
1447 ArrayColumn<Float> specCol(tab, "SPECTRA");
1448 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1449 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
1450 MEpoch::ROScalarColumn timeCol(tab, "TIME");
1451
1452 uInt freqid; freqidCol.get(0, freqid);
1453 Vector<Float> tmpspec; specCol.get(0, tmpspec);
1454 // STAttr.JyPerK has a Vector interface... change sometime.
1455 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
1456 for ( uInt i=0; i<tab.nrow(); ++i) {
1457 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
1458 Float factor = cfac * jyperk;
1459 if ( tokelvin ) factor = Float(1.0) / factor;
1460 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1461 ma *= factor;
1462 specCol.put(i, ma.getArray());
1463 flagCol.put(i, flagsFromMA(ma));
1464 }
1465 ++iter;
1466 }
1467}
1468
1469CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
1470 float tau )
1471{
1472 CountedPtr< Scantable > out = getScantable(in, false);
1473
1474 Table tab = out->table();
1475 ROScalarColumn<Float> elev(tab, "ELEVATION");
1476 ArrayColumn<Float> specCol(tab, "SPECTRA");
1477 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1478 for ( uInt i=0; i<tab.nrow(); ++i) {
1479 Float zdist = Float(C::pi_2) - elev(i);
1480 Float factor = exp(tau/cos(zdist));
1481 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1482 ma *= factor;
1483 specCol.put(i, ma.getArray());
1484 flagCol.put(i, flagsFromMA(ma));
1485 }
1486 return out;
1487}
1488
1489CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
1490 const std::string& kernel,
1491 float width )
1492{
1493 CountedPtr< Scantable > out = getScantable(in, false);
1494 Table& table = out->table();
1495 ArrayColumn<Float> specCol(table, "SPECTRA");
1496 ArrayColumn<uChar> flagCol(table, "FLAGTRA");
1497 Vector<Float> spec;
1498 Vector<uChar> flag;
1499 for ( uInt i=0; i<table.nrow(); ++i) {
1500 specCol.get(i, spec);
1501 flagCol.get(i, flag);
1502 Vector<Bool> mask(flag.nelements());
1503 convertArray(mask, flag);
1504 Vector<Float> specout;
1505 Vector<Bool> maskout;
1506 if ( kernel == "hanning" ) {
1507 mathutil::hanning(specout, maskout, spec , !mask);
1508 convertArray(flag, !maskout);
1509 } else if ( kernel == "rmedian" ) {
1510 mathutil::runningMedian(specout, maskout, spec , mask, width);
1511 convertArray(flag, maskout);
1512 }
1513 flagCol.put(i, flag);
1514 specCol.put(i, specout);
1515 }
1516 return out;
1517}
1518
1519CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
1520 const std::string& kernel, float width )
1521{
1522 if (kernel == "rmedian" || kernel == "hanning") {
1523 return smoothOther(in, kernel, width);
1524 }
1525 CountedPtr< Scantable > out = getScantable(in, false);
1526 Table& table = out->table();
1527 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
1528 // same IFNO should have same no of channels
1529 // this saves overhead
1530 TableIterator iter(table, "IFNO");
1531 while (!iter.pastEnd()) {
1532 Table tab = iter.table();
1533 ArrayColumn<Float> specCol(tab, "SPECTRA");
1534 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1535 Vector<Float> tmpspec; specCol.get(0, tmpspec);
1536 uInt nchan = tmpspec.nelements();
1537 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
1538 Convolver<Float> conv(kvec, IPosition(1,nchan));
1539 Vector<Float> spec;
1540 Vector<uChar> flag;
1541 for ( uInt i=0; i<tab.nrow(); ++i) {
1542 specCol.get(i, spec);
1543 flagCol.get(i, flag);
1544 Vector<Bool> mask(flag.nelements());
1545 convertArray(mask, flag);
1546 Vector<Float> specout;
1547 mathutil::replaceMaskByZero(specout, mask);
1548 conv.linearConv(specout, spec);
1549 specCol.put(i, specout);
1550 }
1551 ++iter;
1552 }
1553 return out;
1554}
1555
1556CountedPtr< Scantable >
1557 STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
1558{
1559 if ( in.size() < 2 ) {
1560 throw(AipsError("Need at least two scantables to perform a merge."));
1561 }
1562 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
1563 bool insitu = insitu_;
1564 setInsitu(false);
1565 CountedPtr< Scantable > out = getScantable(*it, false);
1566 setInsitu(insitu);
1567 Table& tout = out->table();
1568 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
1569 ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
1570 // Renumber SCANNO to be 0-based
1571 Vector<uInt> scannos = scannocol.getColumn();
1572 uInt offset = min(scannos);
1573 scannos -= offset;
1574 scannocol.putColumn(scannos);
1575 uInt newscanno = max(scannos)+1;
1576 ++it;
1577 while ( it != in.end() ){
1578 if ( ! (*it)->conformant(*out) ) {
1579 // log message: "ignoring scantable i, as it isn't
1580 // conformant with the other(s)"
1581 cerr << "oh oh" << endl;
1582 ++it;
1583 continue;
1584 }
1585 out->appendToHistoryTable((*it)->history());
1586 const Table& tab = (*it)->table();
1587 TableIterator scanit(tab, "SCANNO");
1588 while (!scanit.pastEnd()) {
1589 TableIterator freqit(scanit.table(), "FREQ_ID");
1590 while ( !freqit.pastEnd() ) {
1591 Table thetab = freqit.table();
1592 uInt nrow = tout.nrow();
1593 //tout.addRow(thetab.nrow());
1594 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
1595 ROTableRow row(thetab);
1596 for ( uInt i=0; i<thetab.nrow(); ++i) {
1597 uInt k = nrow+i;
1598 scannocol.put(k, newscanno);
1599 const TableRecord& rec = row.get(i);
1600 Double rv,rp,inc;
1601 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
1602 uInt id;
1603 id = out->frequencies().addEntry(rp, rv, inc);
1604 freqidcol.put(k,id);
1605 //String name,fname;Double rf;
1606 Vector<String> name,fname;Vector<Double> rf;
1607 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
1608 id = out->molecules().addEntry(rf, name, fname);
1609 molidcol.put(k, id);
1610 Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
1611 (*it)->focus().getEntry(fax, ftan, frot, fhand,
1612 fmount,fuser, fxy, fxyp,
1613 rec.asuInt("FOCUS_ID"));
1614 id = out->focus().addEntry(fax, ftan, frot, fhand,
1615 fmount,fuser, fxy, fxyp);
1616 focusidcol.put(k, id);
1617 }
1618 ++freqit;
1619 }
1620 ++newscanno;
1621 ++scanit;
1622 }
1623 ++it;
1624 }
1625 return out;
1626}
1627
1628CountedPtr< Scantable >
1629 STMath::invertPhase( const CountedPtr < Scantable >& in )
1630{
1631 return applyToPol(in, &STPol::invertPhase, Float(0.0));
1632}
1633
1634CountedPtr< Scantable >
1635 STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
1636{
1637 return applyToPol(in, &STPol::rotatePhase, Float(phase));
1638}
1639
1640CountedPtr< Scantable >
1641 STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
1642{
1643 return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
1644}
1645
1646CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
1647 STPol::polOperation fptr,
1648 Float phase )
1649{
1650 CountedPtr< Scantable > out = getScantable(in, false);
1651 Table& tout = out->table();
1652 Block<String> cols(4);
1653 cols[0] = String("SCANNO");
1654 cols[1] = String("BEAMNO");
1655 cols[2] = String("IFNO");
1656 cols[3] = String("CYCLENO");
1657 TableIterator iter(tout, cols);
1658 CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
1659 out->getPolType() );
1660 while (!iter.pastEnd()) {
1661 Table t = iter.table();
1662 ArrayColumn<Float> speccol(t, "SPECTRA");
1663 ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
1664 ScalarColumn<Float> parancol(t, "PARANGLE");
1665 Matrix<Float> pols(speccol.getColumn());
1666 try {
1667 stpol->setSpectra(pols);
1668 Float fang,fhand,parang;
1669 fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
1670 fhand = in->focusTable_.getFeedHand(focidcol(0));
1671 parang = parancol(0);
1672 /// @todo re-enable this
1673 // disable total feed angle to support paralactifying Caswell style
1674 stpol->setPhaseCorrections(parang, -parang, fhand);
1675 // use a member function pointer in STPol. This only works on
1676 // the STPol pointer itself, not the Counted Pointer so
1677 // derefernce it.
1678 (&(*(stpol))->*fptr)(phase);
1679 speccol.putColumn(stpol->getSpectra());
1680 } catch (AipsError& e) {
1681 //delete stpol;stpol=0;
1682 throw(e);
1683 }
1684 ++iter;
1685 }
1686 //delete stpol;stpol=0;
1687 return out;
1688}
1689
1690CountedPtr< Scantable >
1691 STMath::swapPolarisations( const CountedPtr< Scantable > & in )
1692{
1693 CountedPtr< Scantable > out = getScantable(in, false);
1694 Table& tout = out->table();
1695 Table t0 = tout(tout.col("POLNO") == 0);
1696 Table t1 = tout(tout.col("POLNO") == 1);
1697 if ( t0.nrow() != t1.nrow() )
1698 throw(AipsError("Inconsistent number of polarisations"));
1699 ArrayColumn<Float> speccol0(t0, "SPECTRA");
1700 ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
1701 ArrayColumn<Float> speccol1(t1, "SPECTRA");
1702 ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
1703 Matrix<Float> s0 = speccol0.getColumn();
1704 Matrix<uChar> f0 = flagcol0.getColumn();
1705 speccol0.putColumn(speccol1.getColumn());
1706 flagcol0.putColumn(flagcol1.getColumn());
1707 speccol1.putColumn(s0);
1708 flagcol1.putColumn(f0);
1709 return out;
1710}
1711
1712CountedPtr< Scantable >
1713 STMath::averagePolarisations( const CountedPtr< Scantable > & in,
1714 const std::vector<bool>& mask,
1715 const std::string& weight )
1716{
1717 if (in->npol() < 2 )
1718 throw(AipsError("averagePolarisations can only be applied to two or more"
1719 "polarisations"));
1720 bool insitu = insitu_;
1721 setInsitu(false);
1722 CountedPtr< Scantable > pols = getScantable(in, true);
1723 setInsitu(insitu);
1724 Table& tout = pols->table();
1725 std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
1726 Table tab = tableCommand(taql, in->table());
1727 if (tab.nrow() == 0 )
1728 throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
1729 TableCopy::copyRows(tout, tab);
1730 TableVector<uInt> vec(tout, "POLNO");
1731 vec = 0;
1732 pols->table_.rwKeywordSet().define("nPol", Int(1));
1733 //pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
1734 pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
1735 std::vector<CountedPtr<Scantable> > vpols;
1736 vpols.push_back(pols);
1737 CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
1738 return out;
1739}
1740
1741CountedPtr< Scantable >
1742 STMath::averageBeams( const CountedPtr< Scantable > & in,
1743 const std::vector<bool>& mask,
1744 const std::string& weight )
1745{
1746 bool insitu = insitu_;
1747 setInsitu(false);
1748 CountedPtr< Scantable > beams = getScantable(in, false);
1749 setInsitu(insitu);
1750 Table& tout = beams->table();
1751 // give all rows the same BEAMNO
1752 TableVector<uInt> vec(tout, "BEAMNO");
1753 vec = 0;
1754 beams->table_.rwKeywordSet().define("nBeam", Int(1));
1755 std::vector<CountedPtr<Scantable> > vbeams;
1756 vbeams.push_back(beams);
1757 CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
1758 return out;
1759}
1760
1761
1762CountedPtr< Scantable >
1763 asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
1764 const std::string & refTime,
1765 const std::string & method)
1766{
1767 // clone as this is not working insitu
1768 bool insitu = insitu_;
1769 setInsitu(false);
1770 CountedPtr< Scantable > out = getScantable(in, false);
1771 setInsitu(insitu);
1772 Table& tout = out->table();
1773 // Get reference Epoch to time of first row or given String
1774 Unit DAY(String("d"));
1775 MEpoch::Ref epochRef(in->getTimeReference());
1776 MEpoch refEpoch;
1777 if (refTime.length()>0) {
1778 Quantum<Double> qt;
1779 if (MVTime::read(qt,refTime)) {
1780 MVEpoch mv(qt);
1781 refEpoch = MEpoch(mv, epochRef);
1782 } else {
1783 throw(AipsError("Invalid format for Epoch string"));
1784 }
1785 } else {
1786 refEpoch = in->timeCol_(0);
1787 }
1788 MPosition refPos = in->getAntennaPosition();
1789
1790 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1791 // test if user frame is different to base frame
1792 if ( in->frequencies().getFrameString(true)
1793 == in->frequencies().getFrameString(false) ) {
1794 throw(AipsError("Can't convert as no output frame has been set"
1795 " (use set_freqframe) or it is aligned already."));
1796 }
1797 MFrequency::Types system = in->frequencies().getFrame();
1798 MVTime mvt(refEpoch.getValue());
1799 String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
1800 ostringstream oss;
1801 oss << "Aligned at reference Epoch " << epochout
1802 << " in frame " << MFrequency::showType(system);
1803 pushLog(String(oss));
1804 // set up the iterator
1805 Block<String> cols(4);
1806 // select by constant direction
1807 cols[0] = String("SRCNAME");
1808 cols[1] = String("BEAMNO");
1809 // select by IF ( no of channels varies over this )
1810 cols[2] = String("IFNO");
1811 // select by restfrequency
1812 cols[3] = String("MOLECULE_ID");
1813 TableIterator iter(tout, cols);
1814 while ( !iter.pastEnd() ) {
1815 Table t = iter.table();
1816 MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1817 TableIterator fiter(t, "FREQ_ID");
1818 // determine nchan from the first row. This should work as
1819 // we are iterating over BEAMNO and IFNO // we should have constant direction
1820
1821 ROArrayColumn<Float> sCol(t, "SPECTRA");
1822 MDirection direction = dirCol(0);
1823 uInt nchan = sCol(0).nelements();
1824 while ( !fiter.pastEnd() ) {
1825 Table ftab = fiter.table();
1826 ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
1827 // get the SpectralCoordinate for the freqid, which we are iterating over
1828 SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1829 FrequencyAligner<Float> fa( sC, nchan, refEpoch,
1830 direction, refPos, system );
1831 // realign the SpectralCoordinate and put into the output Scantable
1832 Vector<String> units(1);
1833 units = String("Hz");
1834 Bool linear=True;
1835 SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1836 sc2.setWorldAxisUnits(units);
1837 uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
1838 sc2.referenceValue()[0],
1839 sc2.increment()[0]);
1840 TableVector<uInt> tvec(ftab, "FREQ_ID");
1841 tvec = id;
1842 // create the "global" abcissa for alignment with same FREQ_ID
1843 Vector<Double> abc(nchan);
1844 Double w;
1845 for (uInt i=0; i<nchan; i++) {
1846 sC.toWorld(w,Double(i));
1847 abc[i] = w;
1848 }
1849 // cache abcissa for same time stamps, so iterate over those
1850 TableIterator timeiter(ftab, "TIME");
1851 while ( !timeiter.pastEnd() ) {
1852 Table tab = timeiter.table();
1853 ArrayColumn<Float> specCol(tab, "SPECTRA");
1854 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1855 MEpoch::ROScalarColumn timeCol(tab, "TIME");
1856 // use align abcissa cache after the first row
1857 bool first = true;
1858 // these rows should be just be POLNO
1859 for (int i=0; i<int(tab.nrow()); ++i) {
1860 // input values
1861 Vector<uChar> flag = flagCol(i);
1862 Vector<Bool> mask(flag.shape());
1863 Vector<Float> specOut, spec;
1864 spec = specCol(i);
1865 Vector<Bool> maskOut;Vector<uChar> flagOut;
1866 convertArray(mask, flag);
1867 // alignment
1868 Bool ok = fa.align(specOut, maskOut, abc, spec,
1869 mask, timeCol(i), !first,
1870 interp, False);
1871 // back into scantable
1872 flagOut.resize(maskOut.nelements());
1873 convertArray(flagOut, maskOut);
1874 flagCol.put(i, flagOut);
1875 specCol.put(i, specOut);
1876 // start abcissa caching
1877 first = false;
1878 }
1879 // next timestamp
1880 ++timeiter;
1881 }
1882 // next FREQ_ID
1883 ++fiter;
1884 }
1885 // next aligner
1886 ++iter;
1887 }
1888 // set this afterwards to ensure we are doing insitu correctly.
1889 out->frequencies().setFrame(system, true);
1890 return out;
1891}
1892
1893CountedPtr<Scantable>
1894 asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
1895 const std::string & newtype )
1896{
1897 if (in->npol() != 2 && in->npol() != 4)
1898 throw(AipsError("Can only convert two or four polarisations."));
1899 if ( in->getPolType() == newtype )
1900 throw(AipsError("No need to convert."));
1901 if ( ! in->selector_.empty() )
1902 throw(AipsError("Can only convert whole scantable. Unset the selection."));
1903 bool insitu = insitu_;
1904 setInsitu(false);
1905 CountedPtr< Scantable > out = getScantable(in, true);
1906 setInsitu(insitu);
1907 Table& tout = out->table();
1908 tout.rwKeywordSet().define("POLTYPE", String(newtype));
1909
1910 Block<String> cols(4);
1911 cols[0] = "SCANNO";
1912 cols[1] = "CYCLENO";
1913 cols[2] = "BEAMNO";
1914 cols[3] = "IFNO";
1915 TableIterator it(in->originalTable_, cols);
1916 String basetype = in->getPolType();
1917 STPol* stpol = STPol::getPolClass(in->factories_, basetype);
1918 try {
1919 while ( !it.pastEnd() ) {
1920 Table tab = it.table();
1921 uInt row = tab.rowNumbers()[0];
1922 stpol->setSpectra(in->getPolMatrix(row));
1923 Float fang,fhand,parang;
1924 fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
1925 fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
1926 parang = in->paraCol_(row);
1927 /// @todo re-enable this
1928 // disable total feed angle to support paralactifying Caswell style
1929 stpol->setPhaseCorrections(parang, -parang, fhand);
1930 Int npolout = 0;
1931 for (uInt i=0; i<tab.nrow(); ++i) {
1932 Vector<Float> outvec = stpol->getSpectrum(i, newtype);
1933 if ( outvec.nelements() > 0 ) {
1934 tout.addRow();
1935 TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
1936 ArrayColumn<Float> sCol(tout,"SPECTRA");
1937 ScalarColumn<uInt> pCol(tout,"POLNO");
1938 sCol.put(tout.nrow()-1 ,outvec);
1939 pCol.put(tout.nrow()-1 ,uInt(npolout));
1940 npolout++;
1941 }
1942 }
1943 tout.rwKeywordSet().define("nPol", npolout);
1944 ++it;
1945 }
1946 } catch (AipsError& e) {
1947 delete stpol;
1948 throw(e);
1949 }
1950 delete stpol;
1951 return out;
1952}
1953
1954CountedPtr< Scantable >
1955 asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
1956 const std::string & scantype )
1957{
1958 bool insitu = insitu_;
1959 setInsitu(false);
1960 CountedPtr< Scantable > out = getScantable(in, true);
1961 setInsitu(insitu);
1962 Table& tout = out->table();
1963 std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
1964 if (scantype == "on") {
1965 taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
1966 }
1967 Table tab = tableCommand(taql, in->table());
1968 TableCopy::copyRows(tout, tab);
1969 if (scantype == "on") {
1970 // re-index SCANNO to 0
1971 TableVector<uInt> vec(tout, "SCANNO");
1972 vec = 0;
1973 }
1974 return out;
1975}
1976
1977CountedPtr< Scantable >
1978 asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
1979 double frequency, double width )
1980{
1981 CountedPtr< Scantable > out = getScantable(in, false);
1982 Table& tout = out->table();
1983 TableIterator iter(tout, "FREQ_ID");
1984 FFTServer<Float,Complex> ffts;
1985 while ( !iter.pastEnd() ) {
1986 Table tab = iter.table();
1987 Double rp,rv,inc;
1988 ROTableRow row(tab);
1989 const TableRecord& rec = row.get(0);
1990 uInt freqid = rec.asuInt("FREQ_ID");
1991 out->frequencies().getEntry(rp, rv, inc, freqid);
1992 ArrayColumn<Float> specCol(tab, "SPECTRA");
1993 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1994 for (int i=0; i<int(tab.nrow()); ++i) {
1995 Vector<Float> spec = specCol(i);
1996 Vector<uChar> flag = flagCol(i);
1997 Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
1998 Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
1999 for (int k=0; k < flag.nelements(); ++k ) {
2000 if (flag[k] > 0) {
2001 spec[k] = 0.0;
2002 }
2003 }
2004 Vector<Complex> lags;
2005 ffts.fft0(lags, spec);
2006 Int start = max(0, lag0);
2007 Int end = min(Int(lags.nelements()-1), lag1);
2008 if (start == end) {
2009 lags[start] = Complex(0.0);
2010 } else {
2011 for (int j=start; j <=end ;++j) {
2012 lags[j] = Complex(0.0);
2013 }
2014 }
2015 ffts.fft0(spec, lags);
2016 specCol.put(i, spec);
2017 }
2018 ++iter;
2019 }
2020 return out;
2021}
2022
2023// Averaging spectra with different channel/resolution
2024CountedPtr<Scantable>
2025STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
2026 const bool& compel,
2027 const std::vector<bool>& mask,
2028 const std::string& weight,
2029 const std::string& avmode )
2030 throw ( casa::AipsError )
2031{
2032 if ( avmode == "SCAN" && in.size() != 1 )
2033 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
2034 "Use merge first."));
2035
2036 CountedPtr<Scantable> out ; // processed result
2037 if ( compel ) {
2038 std::vector< CountedPtr<Scantable> > newin ; // input for average process
2039 uInt insize = in.size() ; // number of scantables
2040
2041 // TEST: do normal average in each table before IF grouping
2042 cout << "Do preliminary averaging" << endl ;
2043 vector< CountedPtr<Scantable> > tmpin( insize ) ;
2044 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2045 vector< CountedPtr<Scantable> > v( 1, in[itable] ) ;
2046 tmpin[itable] = average( v, mask, weight, avmode ) ;
2047 }
2048
2049 // warning
2050 cout << "Average spectra with different spectral resolution" << endl ;
2051 cout << endl ;
2052
2053 // temporarily set coordinfo
2054 vector<string> oldinfo( insize ) ;
2055 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2056 vector<string> coordinfo = in[itable]->getCoordInfo() ;
2057 oldinfo[itable] = coordinfo[0] ;
2058 coordinfo[0] = "Hz" ;
2059 tmpin[itable]->setCoordInfo( coordinfo ) ;
2060 }
2061
2062 // columns
2063 ScalarColumn<uInt> freqIDCol ;
2064 ScalarColumn<uInt> ifnoCol ;
2065 ScalarColumn<uInt> scannoCol ;
2066
2067
2068 // check IF frequency coverage
2069 //cout << "Check IF settings in each table" << endl ;
2070 vector< vector<uInt> > freqid( insize );
2071 vector< vector<double> > iffreq( insize ) ;
2072 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2073 uInt rows = tmpin[itable]->nrow() ;
2074 uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
2075 for ( uInt irow = 0 ; irow < rows ; irow++ ) {
2076 if ( freqid[itable].size() == freqnrows ) {
2077 break ;
2078 }
2079 else {
2080 freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
2081 ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
2082 uInt id = freqIDCol( irow ) ;
2083 if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) {
2084 //cout << "itable = " << itable << ": IF " << id << " is included in the list" << endl ;
2085 vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
2086 freqid[itable].push_back( id ) ;
2087 iffreq[itable].push_back( abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
2088 iffreq[itable].push_back( abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
2089 }
2090 }
2091 }
2092 }
2093
2094 // debug
2095 //cout << "IF settings summary:" << endl ;
2096 //for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
2097 //cout << " Table" << i << endl ;
2098 //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
2099 //cout << " id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
2100 //}
2101 //}
2102 //cout << endl ;
2103
2104 // IF grouping based on their frequency coverage
2105 //cout << "IF grouping based on their frequency coverage" << endl ;
2106
2107 // IF group
2108 // ifgrp[numgrp][nummember*2]
2109 // ifgrp = [table00, freqrow00, table01, freqrow01, ...],
2110 // [table11, freqrow11, table11, freqrow11, ...],
2111 // ...
2112 // ifgfreq[numgrp*2]
2113 // ifgfreq = [min0, max0, min1, max1, ...]
2114 vector< vector<uInt> > ifgrp ;
2115 vector<double> ifgfreq ;
2116
2117 // parameter for IF grouping
2118 // groupmode = OR retrieve all region
2119 // AND only retrieve overlaped region
2120 //string groupmode = "AND" ;
2121 string groupmode = "OR" ;
2122 uInt sizecr = 0 ;
2123 if ( groupmode == "AND" )
2124 sizecr = 2 ;
2125 else if ( groupmode == "OR" )
2126 sizecr = 0 ;
2127
2128 vector<double> sortedfreq ;
2129 for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
2130 for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
2131 if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
2132 sortedfreq.push_back( iffreq[i][j] ) ;
2133 }
2134 }
2135 sort( sortedfreq.begin(), sortedfreq.end() ) ;
2136 for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) {
2137 ifgfreq.push_back( *i ) ;
2138 ifgfreq.push_back( *(i+1) ) ;
2139 }
2140 ifgrp.resize( ifgfreq.size()/2 ) ;
2141 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2142 for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
2143 double range0 = iffreq[itable][2*iif] ;
2144 double range1 = iffreq[itable][2*iif+1] ;
2145 for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) {
2146 double fmin = max( range0, ifgfreq[2*j] ) ;
2147 double fmax = min( range1, ifgfreq[2*j+1] ) ;
2148 if ( fmin < fmax ) {
2149 ifgrp[j].push_back( itable ) ;
2150 ifgrp[j].push_back( freqid[itable][iif] ) ;
2151 }
2152 }
2153 }
2154 }
2155 vector< vector<uInt> >::iterator fiter = ifgrp.begin() ;
2156 vector<double>::iterator giter = ifgfreq.begin() ;
2157 while( fiter != ifgrp.end() ) {
2158 if ( fiter->size() <= sizecr ) {
2159 fiter = ifgrp.erase( fiter ) ;
2160 giter = ifgfreq.erase( giter ) ;
2161 giter = ifgfreq.erase( giter ) ;
2162 }
2163 else {
2164 fiter++ ;
2165 advance( giter, 2 ) ;
2166 }
2167 }
2168
2169 // freqgrp[numgrp][nummember]
2170 // freqgrp = [ifgrp00, ifgrp01, ifgrp02, ...],
2171 // [ifgrp10, ifgrp11, ifgrp12, ...],
2172 // ...
2173 // freqrange[numgrp*2]
2174 // freqrange = [min0, max0, min1, max1, ...]
2175 vector< vector<uInt> > freqgrp ;
2176 double freqrange = 0.0 ;
2177 uInt grpnum = 0 ;
2178 for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
2179 // Assumed that ifgfreq was sorted
2180 if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) {
2181 freqgrp[grpnum-1].push_back( i ) ;
2182 }
2183 else {
2184 vector<uInt> grp0( 1, i ) ;
2185 freqgrp.push_back( grp0 ) ;
2186 grpnum++ ;
2187 }
2188 freqrange = ifgfreq[2*i+1] ;
2189 }
2190
2191
2192 // print IF groups
2193 cout << "IF Group summary: " << endl ;
2194 cout << " GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
2195 for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
2196 cout << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
2197 for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
2198 cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
2199 }
2200 cout << endl ;
2201 }
2202 cout << endl ;
2203
2204 // print frequency group
2205 cout << "Frequency Group summary: " << endl ;
2206 cout << " GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ;
2207 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
2208 cout << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ;
2209 for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) {
2210 cout << freqgrp[i][j] << " " ;
2211 }
2212 cout << endl ;
2213 }
2214 cout << endl ;
2215
2216 // membership check
2217 // groups[numtable][numIF][nummembership]
2218 vector< vector< vector<uInt> > > groups( insize ) ;
2219 for ( uInt i = 0 ; i < insize ; i++ ) {
2220 groups[i].resize( freqid[i].size() ) ;
2221 }
2222 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
2223 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2224 uInt tableid = ifgrp[igrp][2*imem] ;
2225 vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ;
2226 if ( iter != freqid[tableid].end() ) {
2227 uInt rowid = distance( freqid[tableid].begin(), iter ) ;
2228 groups[tableid][rowid].push_back( igrp ) ;
2229 }
2230 }
2231 }
2232
2233 // print membership
2234 //for ( uInt i = 0 ; i < insize ; i++ ) {
2235 //cout << "Table " << i << endl ;
2236 //for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
2237 //cout << " FREQ_ID " << setw( 2 ) << freqid[i][j] << ": " ;
2238 //for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) {
2239 //cout << setw( 2 ) << groups[i][j][k] << " " ;
2240 //}
2241 //cout << endl ;
2242 //}
2243 //}
2244
2245 // set back coordinfo
2246 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2247 vector<string> coordinfo = tmpin[itable]->getCoordInfo() ;
2248 coordinfo[0] = oldinfo[itable] ;
2249 tmpin[itable]->setCoordInfo( coordinfo ) ;
2250 }
2251
2252 // Create additional table if needed
2253 bool oldInsitu = insitu_ ;
2254 setInsitu( false ) ;
2255 vector< vector<uInt> > addrow( insize ) ;
2256 vector<uInt> addtable( insize, 0 ) ;
2257 vector<uInt> newtableids( insize ) ;
2258 vector<uInt> newifids( insize, 0 ) ;
2259 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2260 //cout << "Table " << setw(2) << itable << ": " ;
2261 for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
2262 addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ;
2263 //cout << addrow[itable][ifrow] << " " ;
2264 }
2265 addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ;
2266 //cout << "(" << addtable[itable] << ")" << endl ;
2267 }
2268 newin.resize( insize ) ;
2269 copy( tmpin.begin(), tmpin.end(), newin.begin() ) ;
2270 for ( uInt i = 0 ; i < insize ; i++ ) {
2271 newtableids[i] = i ;
2272 }
2273 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2274 for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
2275 CountedPtr<Scantable> add = getScantable( newin[itable], false ) ;
2276 vector<int> freqidlist ;
2277 for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) {
2278 if ( groups[itable][i].size() > iadd + 1 ) {
2279 freqidlist.push_back( freqid[itable][i] ) ;
2280 }
2281 }
2282 stringstream taqlstream ;
2283 taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ;
2284 for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) {
2285 taqlstream << i ;
2286 if ( i < freqidlist.size() - 1 )
2287 taqlstream << "," ;
2288 else
2289 taqlstream << "]" ;
2290 }
2291 string taql = taqlstream.str() ;
2292 //cout << "taql = " << taql << endl ;
2293 STSelector selector = STSelector() ;
2294 selector.setTaQL( taql ) ;
2295 add->setSelection( selector ) ;
2296 newin.push_back( add ) ;
2297 newtableids.push_back( itable ) ;
2298 newifids.push_back( iadd + 1 ) ;
2299 }
2300 }
2301
2302 // udpate ifgrp
2303 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2304 for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
2305 for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
2306 if ( groups[itable][ifrow].size() > iadd + 1 ) {
2307 uInt igrp = groups[itable][ifrow][iadd+1] ;
2308 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2309 if ( ifgrp[igrp][2*imem] == newtableids[iadd+insize] && ifgrp[igrp][2*imem+1] == freqid[newtableids[iadd+insize]][ifrow] ) {
2310 ifgrp[igrp][2*imem] = insize + iadd ;
2311 }
2312 }
2313 }
2314 }
2315 }
2316 }
2317
2318 // print IF groups again for debug
2319 //cout << "IF Group summary: " << endl ;
2320 //cout << " GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
2321 //for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
2322 //cout << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
2323 //for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
2324 //cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
2325 //}
2326 //cout << endl ;
2327 //}
2328 //cout << endl ;
2329
2330 // reset SCANNO and IF: IF number is reset by the result of sortation
2331 cout << "All scan number is set to 0" << endl ;
2332 //cout << "All IF number is set to IF group index" << endl ;
2333 cout << endl ;
2334 insize = newin.size() ;
2335 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2336 uInt rows = newin[itable]->nrow() ;
2337 Table &tmpt = newin[itable]->table() ;
2338 freqIDCol.attach( tmpt, "FREQ_ID" ) ;
2339 scannoCol.attach( tmpt, "SCANNO" ) ;
2340 ifnoCol.attach( tmpt, "IFNO" ) ;
2341 for ( uInt irow=0 ; irow < rows ; irow++ ) {
2342 scannoCol.put( irow, 0 ) ;
2343 uInt freqID = freqIDCol( irow ) ;
2344 vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ;
2345 if ( iter != freqid[newtableids[itable]].end() ) {
2346 uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ;
2347 ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ;
2348 }
2349 else {
2350 throw(AipsError("IF grouping was wrong in additional tables.")) ;
2351 }
2352 }
2353 }
2354 oldinfo.resize( insize ) ;
2355 setInsitu( oldInsitu ) ;
2356
2357 // temporarily set coordinfo
2358 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2359 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
2360 oldinfo[itable] = coordinfo[0] ;
2361 coordinfo[0] = "Hz" ;
2362 newin[itable]->setCoordInfo( coordinfo ) ;
2363 }
2364
2365 // save column values in the vector
2366 vector< vector<uInt> > freqTableIdVec( insize ) ;
2367 vector< vector<uInt> > freqIdVec( insize ) ;
2368 vector< vector<uInt> > ifNoVec( insize ) ;
2369 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2370 ScalarColumn<uInt> freqIDs ;
2371 freqIDs.attach( newin[itable]->frequencies().table(), "ID" ) ;
2372 ifnoCol.attach( newin[itable]->table(), "IFNO" ) ;
2373 freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ;
2374 for ( uInt irow = 0 ; irow < newin[itable]->frequencies().table().nrow() ; irow++ ) {
2375 freqTableIdVec[itable].push_back( freqIDs( irow ) ) ;
2376 }
2377 for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) {
2378 freqIdVec[itable].push_back( freqIDCol( irow ) ) ;
2379 ifNoVec[itable].push_back( ifnoCol( irow ) ) ;
2380 }
2381 }
2382
2383 // reset spectra and flagtra: pick up common part of frequency coverage
2384 //cout << "Pick common frequency range and align resolution" << endl ;
2385 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2386 uInt rows = newin[itable]->nrow() ;
2387 int nminchan = -1 ;
2388 int nmaxchan = -1 ;
2389 vector<uInt> freqIdUpdate ;
2390 for ( uInt irow = 0 ; irow < rows ; irow++ ) {
2391 uInt ifno = ifNoVec[itable][irow] ; // IFNO is reset by group index
2392 double minfreq = ifgfreq[2*ifno] ;
2393 double maxfreq = ifgfreq[2*ifno+1] ;
2394 //cout << "frequency range: [" << minfreq << "," << maxfreq << "]" << endl ;
2395 vector<double> abcissa = newin[itable]->getAbcissa( irow ) ;
2396 int nchan = abcissa.size() ;
2397 double resol = abcissa[1] - abcissa[0] ;
2398 //cout << "abcissa range : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << endl ;
2399 if ( minfreq <= abcissa[0] )
2400 nminchan = 0 ;
2401 else {
2402 //double cfreq = ( minfreq - abcissa[0] ) / resol ;
2403 double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ;
2404 nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ;
2405 }
2406 if ( maxfreq >= abcissa[abcissa.size()-1] )
2407 nmaxchan = abcissa.size() - 1 ;
2408 else {
2409 //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
2410 double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ;
2411 nmaxchan = abcissa.size() - 1 - int(cfreq) - ( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 ) ;
2412 }
2413 //cout << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << endl ;
2414 if ( nmaxchan > nminchan ) {
2415 newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ;
2416 int newchan = nmaxchan - nminchan + 1 ;
2417 if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan )
2418 freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
2419 }
2420 else {
2421 throw(AipsError("Failed to pick up common part of frequency range.")) ;
2422 }
2423 }
2424 for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
2425 uInt freqId = freqIdUpdate[i] ;
2426 Double refpix ;
2427 Double refval ;
2428 Double increment ;
2429
2430 // update row
2431 newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ;
2432 refval = refval - ( refpix - nminchan ) * increment ;
2433 refpix = 0 ;
2434 newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
2435 }
2436 }
2437
2438
2439 // reset spectra and flagtra: align spectral resolution
2440 //cout << "Align spectral resolution" << endl ;
2441 vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
2442 vector<uInt> gmemid( freqgrp.size(), 0 ) ;
2443 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
2444 double maxdnu = 0.0 ; // maximum (coarsest) frequency resolution
2445 int minchan = INT_MAX ; // minimum channel number
2446 Double refpixref = -1 ; // reference of 'reference pixel'
2447 Double refvalref = -1 ; // reference of 'reference frequency'
2448 Double refinc = -1 ; // reference frequency resolution
2449 uInt refreqid ;
2450 uInt reftable = INT_MAX;
2451 // process only if group member > 1
2452 if ( ifgrp[igrp].size() > 2 ) {
2453 // find minchan and maxdnu in each group
2454 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2455 uInt tableid = ifgrp[igrp][2*imem] ;
2456 uInt rowid = ifgrp[igrp][2*imem+1] ;
2457 vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
2458 if ( iter != freqIdVec[tableid].end() ) {
2459 uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
2460 vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
2461 int nchan = abcissa.size() ;
2462 double dnu = abcissa[1] - abcissa[0] ;
2463 //cout << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << endl ;
2464 if ( nchan < minchan ) {
2465 minchan = nchan ;
2466 maxdnu = dnu ;
2467 newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
2468 refreqid = rowid ;
2469 reftable = tableid ;
2470 }
2471 }
2472 }
2473 // regrid spectra in each group
2474 cout << "GROUP " << igrp << endl ;
2475 cout << " Channel number is adjusted to " << minchan << endl ;
2476 cout << " Corresponding frequency resolution is " << maxdnu << "Hz" << endl ;
2477 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2478 uInt tableid = ifgrp[igrp][2*imem] ;
2479 uInt rowid = ifgrp[igrp][2*imem+1] ;
2480 freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ;
2481 //cout << "tableid = " << tableid << " rowid = " << rowid << ": " << endl ;
2482 //cout << " regridChannel applied to " ;
2483 if ( tableid != reftable )
2484 refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
2485 for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
2486 uInt tfreqid = freqIdVec[tableid][irow] ;
2487 if ( tfreqid == rowid ) {
2488 //cout << irow << " " ;
2489 newin[tableid]->regridChannel( minchan, maxdnu, irow ) ;
2490 freqIDCol.put( irow, refreqid ) ;
2491 freqIdVec[tableid][irow] = refreqid ;
2492 }
2493 }
2494 //cout << endl ;
2495 }
2496 }
2497 else {
2498 uInt tableid = ifgrp[igrp][0] ;
2499 uInt rowid = ifgrp[igrp][1] ;
2500 vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
2501 if ( iter != freqIdVec[tableid].end() ) {
2502 uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
2503 vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
2504 minchan = abcissa.size() ;
2505 maxdnu = abcissa[1] - abcissa[0] ;
2506 }
2507 }
2508 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
2509 if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) {
2510 if ( maxdnu > gmaxdnu[i] ) {
2511 gmaxdnu[i] = maxdnu ;
2512 gmemid[i] = igrp ;
2513 }
2514 break ;
2515 }
2516 }
2517 }
2518 cout << endl ;
2519
2520 // set back coordinfo
2521 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2522 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
2523 coordinfo[0] = oldinfo[itable] ;
2524 newin[itable]->setCoordInfo( coordinfo ) ;
2525 }
2526
2527 // accumulate all rows into the first table
2528 // NOTE: assumed in.size() = 1
2529 vector< CountedPtr<Scantable> > tmp( 1 ) ;
2530 if ( newin.size() == 1 )
2531 tmp[0] = newin[0] ;
2532 else
2533 tmp[0] = merge( newin ) ;
2534
2535 //return tmp[0] ;
2536
2537 // average
2538 CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ;
2539
2540 //return tmpout ;
2541
2542 // combine frequency group
2543 cout << "Combine spectra based on frequency grouping" << endl ;
2544 cout << "IFNO is renumbered as frequency group ID (see above)" << endl ;
2545 vector<string> coordinfo = tmpout->getCoordInfo() ;
2546 oldinfo[0] = coordinfo[0] ;
2547 coordinfo[0] = "Hz" ;
2548 tmpout->setCoordInfo( coordinfo ) ;
2549 // create proformas of output table
2550 stringstream taqlstream ;
2551 taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ;
2552 for ( uInt i = 0 ; i < gmemid.size() ; i++ ) {
2553 taqlstream << gmemid[i] ;
2554 if ( i < gmemid.size() - 1 )
2555 taqlstream << "," ;
2556 else
2557 taqlstream << "]" ;
2558 }
2559 string taql = taqlstream.str() ;
2560 //cout << "taql = " << taql << endl ;
2561 STSelector selector = STSelector() ;
2562 selector.setTaQL( taql ) ;
2563 oldInsitu = insitu_ ;
2564 setInsitu( false ) ;
2565 out = getScantable( tmpout, false ) ;
2566 setInsitu( oldInsitu ) ;
2567 out->setSelection( selector ) ;
2568 // regrid rows
2569 ifnoCol.attach( tmpout->table(), "IFNO" ) ;
2570 for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) {
2571 uInt ifno = ifnoCol( irow ) ;
2572 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
2573 if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) {
2574 vector<double> abcissa = tmpout->getAbcissa( irow ) ;
2575 double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
2576 int nchan = (int)( bw / gmaxdnu[igrp] ) ;
2577 tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
2578 break ;
2579 }
2580 }
2581 }
2582 // combine spectra
2583 ArrayColumn<Float> specColOut ;
2584 specColOut.attach( out->table(), "SPECTRA" ) ;
2585 ArrayColumn<uChar> flagColOut ;
2586 flagColOut.attach( out->table(), "FLAGTRA" ) ;
2587 ScalarColumn<uInt> ifnoColOut ;
2588 ifnoColOut.attach( out->table(), "IFNO" ) ;
2589 ScalarColumn<uInt> polnoColOut ;
2590 polnoColOut.attach( out->table(), "POLNO" ) ;
2591 ScalarColumn<uInt> freqidColOut ;
2592 freqidColOut.attach( out->table(), "FREQ_ID" ) ;
2593 Table &tab = tmpout->table() ;
2594 TableIterator iter( tab, "POLNO" ) ;
2595 vector< vector<uInt> > sizes( freqgrp.size() ) ;
2596 while( !iter.pastEnd() ) {
2597 vector< vector<Float> > specout( freqgrp.size() ) ;
2598 vector< vector<uChar> > flagout( freqgrp.size() ) ;
2599 ArrayColumn<Float> specCols ;
2600 specCols.attach( iter.table(), "SPECTRA" ) ;
2601 ArrayColumn<uChar> flagCols ;
2602 flagCols.attach( iter.table(), "FLAGTRA" ) ;
2603 ifnoCol.attach( iter.table(), "IFNO" ) ;
2604 ScalarColumn<uInt> polnos ;
2605 polnos.attach( iter.table(), "POLNO" ) ;
2606 uInt polno = polnos( 0 ) ;
2607 //cout << "POLNO iteration: " << polno << endl ;
2608 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
2609 sizes[igrp].resize( freqgrp[igrp].size() ) ;
2610 for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
2611 for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
2612 uInt ifno = ifnoCol( irow ) ;
2613 if ( ifno == freqgrp[igrp][imem] ) {
2614 Vector<Float> spec = specCols( irow ) ;
2615 Vector<uChar> flag = flagCols( irow ) ;
2616 vector<Float> svec ;
2617 spec.tovector( svec ) ;
2618 vector<uChar> fvec ;
2619 flag.tovector( fvec ) ;
2620 //cout << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << endl ;
2621 specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
2622 flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
2623 //cout << "specout[" << igrp << "].size() = " << specout[igrp].size() << endl ;
2624 sizes[igrp][imem] = spec.nelements() ;
2625 }
2626 }
2627 }
2628 for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
2629 uInt ifout = ifnoColOut( irow ) ;
2630 uInt polout = polnoColOut( irow ) ;
2631 if ( ifout == gmemid[igrp] && polout == polno ) {
2632 // set SPECTRA and FRAGTRA
2633 Vector<Float> newspec( specout[igrp] ) ;
2634 Vector<uChar> newflag( flagout[igrp] ) ;
2635 specColOut.put( irow, newspec ) ;
2636 flagColOut.put( irow, newflag ) ;
2637 // IFNO renumbering
2638 ifnoColOut.put( irow, igrp ) ;
2639 }
2640 }
2641 }
2642 iter++ ;
2643 }
2644 // update FREQUENCIES subtable
2645 vector<bool> updated( freqgrp.size(), false ) ;
2646 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
2647 uInt index = 0 ;
2648 uInt pixShift = 0 ;
2649 while ( freqgrp[igrp][index] != gmemid[igrp] ) {
2650 pixShift += sizes[igrp][index++] ;
2651 }
2652 for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
2653 if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
2654 uInt freqidOut = freqidColOut( irow ) ;
2655 //cout << "freqgrp " << igrp << " freqidOut = " << freqidOut << endl ;
2656 double refpix ;
2657 double refval ;
2658 double increm ;
2659 out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
2660 refpix += pixShift ;
2661 out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
2662 updated[igrp] = true ;
2663 }
2664 }
2665 }
2666
2667 //out = tmpout ;
2668
2669 coordinfo = tmpout->getCoordInfo() ;
2670 coordinfo[0] = oldinfo[0] ;
2671 tmpout->setCoordInfo( coordinfo ) ;
2672 }
2673 else {
2674 // simple average
2675 out = average( in, mask, weight, avmode ) ;
2676 }
2677
2678 return out ;
2679}
Note: See TracBrowser for help on using the repository browser.