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

Last change on this file since 1481 was 1459, checked in by Takeshi Nakazato, 16 years ago

New Development: No

JIRA Issue: Yes CAS-1084

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: Execute sdcal or sdaverage with scanaverage=True.

Put in Release Notes: No

Description: Fixed a bug in the task sdcal and sdaverage

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