source: trunk/src/STMath.cpp@ 1180

Last change on this file since 1180 was 1145, checked in by mar637, 18 years ago

added average_beam which is part of Ticket #45
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 45.2 KB
Line 
1//
2// C++ Implementation: STMath
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13#include <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
34#include <lattices/Lattices/LatticeUtilities.h>
35
36#include <coordinates/Coordinates/SpectralCoordinate.h>
37#include <coordinates/Coordinates/CoordinateSystem.h>
38#include <coordinates/Coordinates/CoordinateUtil.h>
39#include <coordinates/Coordinates/FrequencyAligner.h>
40
41#include <scimath/Mathematics/VectorKernel.h>
42#include <scimath/Mathematics/Convolver.h>
43#include <scimath/Functionals/Polynomial.h>
44
45#include "MathUtils.h"
46#include "RowAccumulator.h"
47#include "STAttr.h"
48#include "STMath.h"
49
50using namespace casa;
51
52using namespace asap;
53
54STMath::STMath(bool insitu) :
55 insitu_(insitu)
56{
57}
58
59
60STMath::~STMath()
61{
62}
63
64CountedPtr<Scantable>
65STMath::average( const std::vector<CountedPtr<Scantable> >& in,
66 const std::vector<bool>& mask,
67 const std::string& weight,
68 const std::string& avmode)
69{
70 if ( avmode == "SCAN" && in.size() != 1 )
71 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
72 "Use merge first."));
73 WeightType wtype = stringToWeight(weight);
74
75 // output
76 // clone as this is non insitu
77 bool insitu = insitu_;
78 setInsitu(false);
79 CountedPtr< Scantable > out = getScantable(in[0], true);
80 setInsitu(insitu);
81 std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
82 ++stit;
83 while ( stit != in.end() ) {
84 out->appendToHistoryTable((*stit)->history());
85 ++stit;
86 }
87
88 Table& tout = out->table();
89
90 /// @todo check if all scantables are conformant
91
92 ArrayColumn<Float> specColOut(tout,"SPECTRA");
93 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
94 ArrayColumn<Float> tsysColOut(tout,"TSYS");
95 ScalarColumn<Double> mjdColOut(tout,"TIME");
96 ScalarColumn<Double> intColOut(tout,"INTERVAL");
97 ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
98 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
99
100 // set up the output table rows. These are based on the structure of the
101 // FIRST scantable in the vector
102 const Table& baset = in[0]->table();
103
104 Block<String> cols(3);
105 cols[0] = String("BEAMNO");
106 cols[1] = String("IFNO");
107 cols[2] = String("POLNO");
108 if ( avmode == "SOURCE" ) {
109 cols.resize(4);
110 cols[3] = String("SRCNAME");
111 }
112 if ( avmode == "SCAN" && in.size() == 1) {
113 cols.resize(4);
114 cols[3] = String("SCANNO");
115 }
116 uInt outrowCount = 0;
117 TableIterator iter(baset, cols);
118 while (!iter.pastEnd()) {
119 Table subt = iter.table();
120 // copy the first row of this selection into the new table
121 tout.addRow();
122 TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
123 // re-index to 0
124 if ( avmode != "SCAN" && avmode != "SOURCE" ) {
125 scanColOut.put(outrowCount, uInt(0));
126 }
127 ++outrowCount;
128 ++iter;
129 }
130 RowAccumulator acc(wtype);
131 Vector<Bool> cmask(mask);
132 acc.setUserMask(cmask);
133 ROTableRow row(tout);
134 ROArrayColumn<Float> specCol, tsysCol;
135 ROArrayColumn<uChar> flagCol;
136 ROScalarColumn<Double> mjdCol, intCol;
137 ROScalarColumn<Int> scanIDCol;
138
139 for (uInt i=0; i < tout.nrow(); ++i) {
140 for ( int j=0; j < int(in.size()); ++j ) {
141 const Table& tin = in[j]->table();
142 const TableRecord& rec = row.get(i);
143 ROScalarColumn<Double> tmp(tin, "TIME");
144 Double td;tmp.get(0,td);
145 Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
146 && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
147 && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
148 Table subt;
149 if ( avmode == "SOURCE") {
150 subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
151 } else if (avmode == "SCAN") {
152 subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
153 } else {
154 subt = basesubt;
155 }
156 specCol.attach(subt,"SPECTRA");
157 flagCol.attach(subt,"FLAGTRA");
158 tsysCol.attach(subt,"TSYS");
159 intCol.attach(subt,"INTERVAL");
160 mjdCol.attach(subt,"TIME");
161 Vector<Float> spec,tsys;
162 Vector<uChar> flag;
163 Double inter,time;
164 for (uInt k = 0; k < subt.nrow(); ++k ) {
165 flagCol.get(k, flag);
166 Vector<Bool> bflag(flag.shape());
167 convertArray(bflag, flag);
168 if ( allEQ(bflag, True) ) {
169 continue;//don't accumulate
170 }
171 specCol.get(k, spec);
172 tsysCol.get(k, tsys);
173 intCol.get(k, inter);
174 mjdCol.get(k, time);
175 // spectrum has to be added last to enable weighting by the other values
176 acc.add(spec, !bflag, tsys, inter, time);
177 }
178 }
179 //write out
180 specColOut.put(i, acc.getSpectrum());
181 const Vector<Bool>& msk = acc.getMask();
182 Vector<uChar> flg(msk.shape());
183 convertArray(flg, !msk);
184 flagColOut.put(i, flg);
185 tsysColOut.put(i, acc.getTsys());
186 intColOut.put(i, acc.getInterval());
187 mjdColOut.put(i, acc.getTime());
188 // we should only have one cycle now -> reset it to be 0
189 // frequency switched data has different CYCLENO for different IFNO
190 // which requires resetting this value
191 cycColOut.put(i, uInt(0));
192 acc.reset();
193 }
194 return out;
195}
196
197CountedPtr< Scantable >
198 STMath::averageChannel( const CountedPtr < Scantable > & in,
199 const std::string & mode,
200 const std::string& avmode )
201{
202 // clone as this is non insitu
203 bool insitu = insitu_;
204 setInsitu(false);
205 CountedPtr< Scantable > out = getScantable(in, true);
206 setInsitu(insitu);
207 Table& tout = out->table();
208 ArrayColumn<Float> specColOut(tout,"SPECTRA");
209 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
210 ArrayColumn<Float> tsysColOut(tout,"TSYS");
211 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
212 Table tmp = in->table().sort("BEAMNO");
213 Block<String> cols(3);
214 cols[0] = String("BEAMNO");
215 cols[1] = String("IFNO");
216 cols[2] = String("POLNO");
217 if ( avmode == "SCAN") {
218 cols.resize(4);
219 cols[3] = String("SCANNO");
220 }
221 uInt outrowCount = 0;
222 uChar userflag = 1 << 7;
223 TableIterator iter(tmp, cols);
224 while (!iter.pastEnd()) {
225 Table subt = iter.table();
226 ROArrayColumn<Float> specCol, tsysCol;
227 ROArrayColumn<uChar> flagCol;
228 specCol.attach(subt,"SPECTRA");
229 flagCol.attach(subt,"FLAGTRA");
230 tsysCol.attach(subt,"TSYS");
231 tout.addRow();
232 TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
233 if ( avmode != "SCAN") {
234 scanColOut.put(outrowCount, uInt(0));
235 }
236 Vector<Float> tmp;
237 specCol.get(0, tmp);
238 uInt nchan = tmp.nelements();
239 // have to do channel by channel here as MaskedArrMath
240 // doesn't have partialMedians
241 Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
242 Vector<Float> outspec(nchan);
243 Vector<uChar> outflag(nchan,0);
244 Vector<Float> outtsys(1);/// @fixme when tsys is channel based
245 for (uInt i=0; i<nchan; ++i) {
246 Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
247 MaskedArray<Float> ma = maskedArray(specs,flags);
248 outspec[i] = median(ma);
249 if ( allEQ(ma.getMask(), False) )
250 outflag[i] = userflag;// flag data
251 }
252 outtsys[0] = median(tsysCol.getColumn());
253 specColOut.put(outrowCount, outspec);
254 flagColOut.put(outrowCount, outflag);
255 tsysColOut.put(outrowCount, outtsys);
256
257 ++outrowCount;
258 ++iter;
259 }
260 return out;
261}
262
263CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
264 bool droprows)
265{
266 if (insitu_) return in;
267 else {
268 // clone
269 Scantable* tabp = new Scantable(*in, Bool(droprows));
270 return CountedPtr<Scantable>(tabp);
271 }
272}
273
274CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
275 float val,
276 const std::string& mode,
277 bool tsys )
278{
279 // modes are "ADD" and "MUL"
280 CountedPtr< Scantable > out = getScantable(in, false);
281 Table& tab = out->table();
282 ArrayColumn<Float> specCol(tab,"SPECTRA");
283 ArrayColumn<Float> tsysCol(tab,"TSYS");
284 for (uInt i=0; i<tab.nrow(); ++i) {
285 Vector<Float> spec;
286 Vector<Float> ts;
287 specCol.get(i, spec);
288 tsysCol.get(i, ts);
289 if (mode == "MUL") {
290 spec *= val;
291 specCol.put(i, spec);
292 if ( tsys ) {
293 ts *= val;
294 tsysCol.put(i, ts);
295 }
296 } else if ( mode == "ADD" ) {
297 spec += val;
298 specCol.put(i, spec);
299 if ( tsys ) {
300 ts += val;
301 tsysCol.put(i, ts);
302 }
303 }
304 }
305 return out;
306}
307
308MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
309 const Vector<uChar>& f)
310{
311 Vector<Bool> mask;
312 mask.resize(f.shape());
313 convertArray(mask, f);
314 return MaskedArray<Float>(s,!mask);
315}
316
317Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
318{
319 const Vector<Bool>& m = ma.getMask();
320 Vector<uChar> flags(m.shape());
321 convertArray(flags, !m);
322 return flags;
323}
324
325CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
326 const std::string & mode,
327 bool preserve )
328{
329 /// @todo make other modes available
330 /// modes should be "nearest", "pair"
331 // make this operation non insitu
332 const Table& tin = in->table();
333 Table ons = tin(tin.col("SRCTYPE") == Int(0));
334 Table offs = tin(tin.col("SRCTYPE") == Int(1));
335 if ( offs.nrow() == 0 )
336 throw(AipsError("No 'off' scans present."));
337 // put all "on" scans into output table
338
339 bool insitu = insitu_;
340 setInsitu(false);
341 CountedPtr< Scantable > out = getScantable(in, true);
342 setInsitu(insitu);
343 Table& tout = out->table();
344
345 TableCopy::copyRows(tout, ons);
346 TableRow row(tout);
347 ROScalarColumn<Double> offtimeCol(offs, "TIME");
348
349 ArrayColumn<Float> outspecCol(tout, "SPECTRA");
350 ROArrayColumn<Float> outtsysCol(tout, "TSYS");
351 ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
352 for (uInt i=0; i < tout.nrow(); ++i) {
353 const TableRecord& rec = row.get(i);
354 Double ontime = rec.asDouble("TIME");
355 ROScalarColumn<Double> offtimeCol(offs, "TIME");
356 Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
357 Table sel = offs( abs(offs.col("TIME")-ontime) <= mindeltat
358 && offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
359 && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
360 && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
361
362 TableRow offrow(sel);
363 const TableRecord& offrec = offrow.get(0);//should only be one row
364 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
365 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
366 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
367 /// @fixme this assumes tsys is a scalar not vector
368 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
369 Vector<Float> specon, tsyson;
370 outtsysCol.get(i, tsyson);
371 outspecCol.get(i, specon);
372 Vector<uChar> flagon;
373 outflagCol.get(i, flagon);
374 MaskedArray<Float> mon = maskedArray(specon, flagon);
375 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
376 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
377 if (preserve) {
378 quot -= tsysoffscalar;
379 } else {
380 quot -= tsyson[0];
381 }
382 outspecCol.put(i, quot.getArray());
383 outflagCol.put(i, flagsFromMA(quot));
384 }
385 // renumber scanno
386 TableIterator it(tout, "SCANNO");
387 uInt i = 0;
388 while ( !it.pastEnd() ) {
389 Table t = it.table();
390 TableVector<uInt> vec(t, "SCANNO");
391 vec = i;
392 ++i;
393 ++it;
394 }
395 return out;
396}
397
398
399CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
400 const CountedPtr< Scantable > & off,
401 bool preserve )
402{
403 bool insitu = insitu_;
404 if ( ! on->conformant(*off) ) {
405 throw(AipsError("'on' and 'off' scantables are not conformant."));
406 }
407 setInsitu(false);
408 CountedPtr< Scantable > out = getScantable(on, false);
409 setInsitu(insitu);
410 Table& tout = out->table();
411 const Table& toff = off->table();
412 TableIterator sit(tout, "SCANNO");
413 TableIterator s2it(toff, "SCANNO");
414 while ( !sit.pastEnd() ) {
415 Table ton = sit.table();
416 TableRow row(ton);
417 Table t = s2it.table();
418 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
419 ROArrayColumn<Float> outtsysCol(ton, "TSYS");
420 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
421 for (uInt i=0; i < ton.nrow(); ++i) {
422 const TableRecord& rec = row.get(i);
423 Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
424 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
425 && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
426 if ( offsel.nrow() == 0 )
427 throw AipsError("STMath::quotient: no matching off");
428 TableRow offrow(offsel);
429 const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
430 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
431 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
432 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
433 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
434 Vector<Float> specon, tsyson;
435 outtsysCol.get(i, tsyson);
436 outspecCol.get(i, specon);
437 Vector<uChar> flagon;
438 outflagCol.get(i, flagon);
439 MaskedArray<Float> mon = maskedArray(specon, flagon);
440 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
441 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
442 if (preserve) {
443 quot -= tsysoffscalar;
444 } else {
445 quot -= tsyson[0];
446 }
447 outspecCol.put(i, quot.getArray());
448 outflagCol.put(i, flagsFromMA(quot));
449 }
450 ++sit;
451 ++s2it;
452 // take the first off for each on scan which doesn't have a
453 // matching off scan
454 // non <= noff: matching pairs, non > noff matching pairs then first off
455 if ( s2it.pastEnd() ) s2it.reset();
456 }
457 return out;
458}
459
460
461CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
462{
463 // make copy or reference
464 CountedPtr< Scantable > out = getScantable(in, false);
465 Table& tout = out->table();
466 Block<String> cols(4);
467 cols[0] = String("SCANNO");
468 cols[1] = String("CYCLENO");
469 cols[2] = String("BEAMNO");
470 cols[3] = String("POLNO");
471 TableIterator iter(tout, cols);
472 while (!iter.pastEnd()) {
473 Table subt = iter.table();
474 // this should leave us with two rows for the two IFs....if not ignore
475 if (subt.nrow() != 2 ) {
476 continue;
477 }
478 ArrayColumn<Float> specCol(subt, "SPECTRA");
479 ArrayColumn<Float> tsysCol(subt, "TSYS");
480 ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
481 Vector<Float> onspec,offspec, ontsys, offtsys;
482 Vector<uChar> onflag, offflag;
483 tsysCol.get(0, ontsys); tsysCol.get(1, offtsys);
484 specCol.get(0, onspec); specCol.get(1, offspec);
485 flagCol.get(0, onflag); flagCol.get(1, offflag);
486 MaskedArray<Float> on = maskedArray(onspec, onflag);
487 MaskedArray<Float> off = maskedArray(offspec, offflag);
488 MaskedArray<Float> oncopy = on.copy();
489
490 on /= off; on -= 1.0f;
491 on *= ontsys[0];
492 off /= oncopy; off -= 1.0f;
493 off *= offtsys[0];
494 specCol.put(0, on.getArray());
495 const Vector<Bool>& m0 = on.getMask();
496 Vector<uChar> flags0(m0.shape());
497 convertArray(flags0, !m0);
498 flagCol.put(0, flags0);
499
500 specCol.put(1, off.getArray());
501 const Vector<Bool>& m1 = off.getMask();
502 Vector<uChar> flags1(m1.shape());
503 convertArray(flags1, !m1);
504 flagCol.put(1, flags1);
505 ++iter;
506 }
507
508 return out;
509}
510
511std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
512 const std::vector< bool > & mask,
513 const std::string& which )
514{
515
516 Vector<Bool> m(mask);
517 const Table& tab = in->table();
518 ROArrayColumn<Float> specCol(tab, "SPECTRA");
519 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
520 std::vector<float> out;
521 for (uInt i=0; i < tab.nrow(); ++i ) {
522 Vector<Float> spec; specCol.get(i, spec);
523 Vector<uChar> flag; flagCol.get(i, flag);
524 MaskedArray<Float> ma = maskedArray(spec, flag);
525 float outstat = 0.0;
526 if ( spec.nelements() == m.nelements() ) {
527 outstat = mathutil::statistics(which, ma(m));
528 } else {
529 outstat = mathutil::statistics(which, ma);
530 }
531 out.push_back(outstat);
532 }
533 return out;
534}
535
536CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
537 int width )
538{
539 if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
540 CountedPtr< Scantable > out = getScantable(in, false);
541 Table& tout = out->table();
542 out->frequencies().rescale(width, "BIN");
543 ArrayColumn<Float> specCol(tout, "SPECTRA");
544 ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
545 for (uInt i=0; i < tout.nrow(); ++i ) {
546 MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
547 MaskedArray<Float> maout;
548 LatticeUtilities::bin(maout, main, 0, Int(width));
549 /// @todo implement channel based tsys binning
550 specCol.put(i, maout.getArray());
551 flagCol.put(i, flagsFromMA(maout));
552 // take only the first binned spectrum's length for the deprecated
553 // global header item nChan
554 if (i==0) tout.rwKeywordSet().define(String("nChan"),
555 Int(maout.getArray().nelements()));
556 }
557 return out;
558}
559
560CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
561 const std::string& method,
562 float width )
563//
564// Should add the possibility of width being specified in km/s. This means
565// that for each freqID (SpectralCoordinate) we will need to convert to an
566// average channel width (say at the reference pixel). Then we would need
567// to be careful to make sure each spectrum (of different freqID)
568// is the same length.
569//
570{
571 //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
572 Int interpMethod(stringToIMethod(method));
573
574 CountedPtr< Scantable > out = getScantable(in, false);
575 Table& tout = out->table();
576
577// Resample SpectralCoordinates (one per freqID)
578 out->frequencies().rescale(width, "RESAMPLE");
579 TableIterator iter(tout, "IFNO");
580 TableRow row(tout);
581 while ( !iter.pastEnd() ) {
582 Table tab = iter.table();
583 ArrayColumn<Float> specCol(tab, "SPECTRA");
584 //ArrayColumn<Float> tsysCol(tout, "TSYS");
585 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
586 Vector<Float> spec;
587 Vector<uChar> flag;
588 specCol.get(0,spec); // the number of channels should be constant per IF
589 uInt nChanIn = spec.nelements();
590 Vector<Float> xIn(nChanIn); indgen(xIn);
591 Int fac = Int(nChanIn/width);
592 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
593 uInt k = 0;
594 Float x = 0.0;
595 while (x < Float(nChanIn) ) {
596 xOut(k) = x;
597 k++;
598 x += width;
599 }
600 uInt nChanOut = k;
601 xOut.resize(nChanOut, True);
602 // process all rows for this IFNO
603 Vector<Float> specOut;
604 Vector<Bool> maskOut;
605 Vector<uChar> flagOut;
606 for (uInt i=0; i < tab.nrow(); ++i) {
607 specCol.get(i, spec);
608 flagCol.get(i, flag);
609 Vector<Bool> mask(flag.nelements());
610 convertArray(mask, flag);
611
612 IPosition shapeIn(spec.shape());
613 //sh.nchan = nChanOut;
614 InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
615 xIn, spec, mask,
616 interpMethod, True, True);
617 /// @todo do the same for channel based Tsys
618 flagOut.resize(maskOut.nelements());
619 convertArray(flagOut, maskOut);
620 specCol.put(i, specOut);
621 flagCol.put(i, flagOut);
622 }
623 ++iter;
624 }
625
626 return out;
627}
628
629STMath::imethod STMath::stringToIMethod(const std::string& in)
630{
631 static STMath::imap lookup;
632
633 // initialize the lookup table if necessary
634 if ( lookup.empty() ) {
635 lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
636 lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
637 lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
638 lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
639 }
640
641 STMath::imap::const_iterator iter = lookup.find(in);
642
643 if ( lookup.end() == iter ) {
644 std::string message = in;
645 message += " is not a valid interpolation mode";
646 throw(AipsError(message));
647 }
648 return iter->second;
649}
650
651WeightType STMath::stringToWeight(const std::string& in)
652{
653 static std::map<std::string, WeightType> lookup;
654
655 // initialize the lookup table if necessary
656 if ( lookup.empty() ) {
657 lookup["NONE"] = asap::NONE;
658 lookup["TINT"] = asap::TINT;
659 lookup["TINTSYS"] = asap::TINTSYS;
660 lookup["TSYS"] = asap::TSYS;
661 lookup["VAR"] = asap::VAR;
662 }
663
664 std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
665
666 if ( lookup.end() == iter ) {
667 std::string message = in;
668 message += " is not a valid weighting mode";
669 throw(AipsError(message));
670 }
671 return iter->second;
672}
673
674CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
675 const vector< float > & coeff,
676 const std::string & filename,
677 const std::string& method)
678{
679 // Get elevation data from Scantable and convert to degrees
680 CountedPtr< Scantable > out = getScantable(in, false);
681 Table& tab = out->table();
682 ROScalarColumn<Float> elev(tab, "ELEVATION");
683 Vector<Float> x = elev.getColumn();
684 x *= Float(180 / C::pi); // Degrees
685
686 Vector<Float> coeffs(coeff);
687 const uInt nc = coeffs.nelements();
688 if ( filename.length() > 0 && nc > 0 ) {
689 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
690 }
691
692 // Correct
693 if ( nc > 0 || filename.length() == 0 ) {
694 // Find instrument
695 Bool throwit = True;
696 Instrument inst =
697 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
698 throwit);
699
700 // Set polynomial
701 Polynomial<Float>* ppoly = 0;
702 Vector<Float> coeff;
703 String msg;
704 if ( nc > 0 ) {
705 ppoly = new Polynomial<Float>(nc);
706 coeff = coeffs;
707 msg = String("user");
708 } else {
709 STAttr sdAttr;
710 coeff = sdAttr.gainElevationPoly(inst);
711 ppoly = new Polynomial<Float>(3);
712 msg = String("built in");
713 }
714
715 if ( coeff.nelements() > 0 ) {
716 ppoly->setCoefficients(coeff);
717 } else {
718 delete ppoly;
719 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
720 }
721 ostringstream oss;
722 oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
723 oss << " " << coeff;
724 pushLog(String(oss));
725 const uInt nrow = tab.nrow();
726 Vector<Float> factor(nrow);
727 for ( uInt i=0; i < nrow; ++i ) {
728 factor[i] = 1.0 / (*ppoly)(x[i]);
729 }
730 delete ppoly;
731 scaleByVector(tab, factor, true);
732
733 } else {
734 // Read and correct
735 pushLog("Making correction from ascii Table");
736 scaleFromAsciiTable(tab, filename, method, x, true);
737 }
738 return out;
739}
740
741void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
742 const std::string& method,
743 const Vector<Float>& xout, bool dotsys)
744{
745
746// Read gain-elevation ascii file data into a Table.
747
748 String formatString;
749 Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
750 scaleFromTable(in, tbl, method, xout, dotsys);
751}
752
753void STMath::scaleFromTable(Table& in,
754 const Table& table,
755 const std::string& method,
756 const Vector<Float>& xout, bool dotsys)
757{
758
759 ROScalarColumn<Float> geElCol(table, "ELEVATION");
760 ROScalarColumn<Float> geFacCol(table, "FACTOR");
761 Vector<Float> xin = geElCol.getColumn();
762 Vector<Float> yin = geFacCol.getColumn();
763 Vector<Bool> maskin(xin.nelements(),True);
764
765 // Interpolate (and extrapolate) with desired method
766
767 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
768
769 Vector<Float> yout;
770 Vector<Bool> maskout;
771 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
772 xin, yin, maskin, interp,
773 True, True);
774
775 scaleByVector(in, Float(1.0)/yout, dotsys);
776}
777
778void STMath::scaleByVector( Table& in,
779 const Vector< Float >& factor,
780 bool dotsys )
781{
782 uInt nrow = in.nrow();
783 if ( factor.nelements() != nrow ) {
784 throw(AipsError("factors.nelements() != table.nelements()"));
785 }
786 ArrayColumn<Float> specCol(in, "SPECTRA");
787 ArrayColumn<uChar> flagCol(in, "FLAGTRA");
788 ArrayColumn<Float> tsysCol(in, "TSYS");
789 for (uInt i=0; i < nrow; ++i) {
790 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
791 ma *= factor[i];
792 specCol.put(i, ma.getArray());
793 flagCol.put(i, flagsFromMA(ma));
794 if ( dotsys ) {
795 Vector<Float> tsys = tsysCol(i);
796 tsys *= factor[i];
797 tsysCol.put(i,tsys);
798 }
799 }
800}
801
802CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
803 float d, float etaap,
804 float jyperk )
805{
806 CountedPtr< Scantable > out = getScantable(in, false);
807 Table& tab = in->table();
808 Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
809 Unit K(String("K"));
810 Unit JY(String("Jy"));
811
812 bool tokelvin = true;
813 Double cfac = 1.0;
814
815 if ( fluxUnit == JY ) {
816 pushLog("Converting to K");
817 Quantum<Double> t(1.0,fluxUnit);
818 Quantum<Double> t2 = t.get(JY);
819 cfac = (t2 / t).getValue(); // value to Jy
820
821 tokelvin = true;
822 out->setFluxUnit("K");
823 } else if ( fluxUnit == K ) {
824 pushLog("Converting to Jy");
825 Quantum<Double> t(1.0,fluxUnit);
826 Quantum<Double> t2 = t.get(K);
827 cfac = (t2 / t).getValue(); // value to K
828
829 tokelvin = false;
830 out->setFluxUnit("Jy");
831 } else {
832 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
833 }
834 // Make sure input values are converted to either Jy or K first...
835 Float factor = cfac;
836
837 // Select method
838 if (jyperk > 0.0) {
839 factor *= jyperk;
840 if ( tokelvin ) factor = 1.0 / jyperk;
841 ostringstream oss;
842 oss << "Jy/K = " << jyperk;
843 pushLog(String(oss));
844 Vector<Float> factors(tab.nrow(), factor);
845 scaleByVector(tab,factors, false);
846 } else if ( etaap > 0.0) {
847 Instrument inst =
848 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), True);
849 STAttr sda;
850 if (d < 0) d = sda.diameter(inst);
851 jyperk = STAttr::findJyPerK(etaap, d);
852 ostringstream oss;
853 oss << "Jy/K = " << jyperk;
854 pushLog(String(oss));
855 factor *= jyperk;
856 if ( tokelvin ) {
857 factor = 1.0 / factor;
858 }
859 Vector<Float> factors(tab.nrow(), factor);
860 scaleByVector(tab, factors, False);
861 } else {
862
863 // OK now we must deal with automatic look up of values.
864 // We must also deal with the fact that the factors need
865 // to be computed per IF and may be different and may
866 // change per integration.
867
868 pushLog("Looking up conversion factors");
869 convertBrightnessUnits(out, tokelvin, cfac);
870 }
871
872 return out;
873}
874
875void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
876 bool tokelvin, float cfac )
877{
878 Table& table = in->table();
879 Instrument inst =
880 STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
881 TableIterator iter(table, "FREQ_ID");
882 STFrequencies stfreqs = in->frequencies();
883 STAttr sdAtt;
884 while (!iter.pastEnd()) {
885 Table tab = iter.table();
886 ArrayColumn<Float> specCol(tab, "SPECTRA");
887 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
888 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
889 MEpoch::ROScalarColumn timeCol(tab, "TIME");
890
891 uInt freqid; freqidCol.get(0, freqid);
892 Vector<Float> tmpspec; specCol.get(0, tmpspec);
893 // STAttr.JyPerK has a Vector interface... change sometime.
894 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
895 for ( uInt i=0; i<tab.nrow(); ++i) {
896 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
897 Float factor = cfac * jyperk;
898 if ( tokelvin ) factor = Float(1.0) / factor;
899 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
900 ma *= factor;
901 specCol.put(i, ma.getArray());
902 flagCol.put(i, flagsFromMA(ma));
903 }
904 ++iter;
905 }
906}
907
908CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
909 float tau )
910{
911 CountedPtr< Scantable > out = getScantable(in, false);
912
913 Table tab = out->table();
914 ROScalarColumn<Float> elev(tab, "ELEVATION");
915 ArrayColumn<Float> specCol(tab, "SPECTRA");
916 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
917 for ( uInt i=0; i<tab.nrow(); ++i) {
918 Float zdist = Float(C::pi_2) - elev(i);
919 Float factor = exp(tau)/cos(zdist);
920 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
921 ma *= factor;
922 specCol.put(i, ma.getArray());
923 flagCol.put(i, flagsFromMA(ma));
924 }
925 return out;
926}
927
928CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
929 const std::string& kernel, float width )
930{
931 CountedPtr< Scantable > out = getScantable(in, false);
932 Table& table = out->table();
933 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
934 // same IFNO should have same no of channels
935 // this saves overhead
936 TableIterator iter(table, "IFNO");
937 while (!iter.pastEnd()) {
938 Table tab = iter.table();
939 ArrayColumn<Float> specCol(tab, "SPECTRA");
940 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
941 Vector<Float> tmpspec; specCol.get(0, tmpspec);
942 uInt nchan = tmpspec.nelements();
943 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
944 Convolver<Float> conv(kvec, IPosition(1,nchan));
945 Vector<Float> spec;
946 Vector<uChar> flag;
947 for ( uInt i=0; i<tab.nrow(); ++i) {
948 specCol.get(i, spec);
949 flagCol.get(i, flag);
950 Vector<Bool> mask(flag.nelements());
951 convertArray(mask, flag);
952 Vector<Float> specout;
953 if ( type == VectorKernel::HANNING ) {
954 Vector<Bool> maskout;
955 mathutil::hanning(specout, maskout, spec , !mask);
956 convertArray(flag, !maskout);
957 flagCol.put(i, flag);
958 specCol.put(i, specout);
959 } else {
960 mathutil::replaceMaskByZero(specout, mask);
961 conv.linearConv(specout, spec);
962 specCol.put(i, specout);
963 }
964 }
965 ++iter;
966 }
967 return out;
968}
969
970CountedPtr< Scantable >
971 STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
972{
973 if ( in.size() < 2 ) {
974 throw(AipsError("Need at least two scantables to perform a merge."));
975 }
976 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
977 bool insitu = insitu_;
978 setInsitu(false);
979 CountedPtr< Scantable > out = getScantable(*it, false);
980 setInsitu(insitu);
981 Table& tout = out->table();
982 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
983 ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
984 // Renumber SCANNO to be 0-based
985 Vector<uInt> scannos = scannocol.getColumn();
986 uInt offset = min(scannos);
987 scannos -= offset;
988 scannocol.putColumn(scannos);
989 uInt newscanno = max(scannos)+1;
990 ++it;
991 while ( it != in.end() ){
992 if ( ! (*it)->conformant(*out) ) {
993 // log message: "ignoring scantable i, as it isn't
994 // conformant with the other(s)"
995 cerr << "oh oh" << endl;
996 ++it;
997 continue;
998 }
999 out->appendToHistoryTable((*it)->history());
1000 const Table& tab = (*it)->table();
1001 TableIterator scanit(tab, "SCANNO");
1002 while (!scanit.pastEnd()) {
1003 TableIterator freqit(scanit.table(), "FREQ_ID");
1004 while ( !freqit.pastEnd() ) {
1005 Table thetab = freqit.table();
1006 uInt nrow = tout.nrow();
1007 //tout.addRow(thetab.nrow());
1008 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
1009 ROTableRow row(thetab);
1010 for ( uInt i=0; i<thetab.nrow(); ++i) {
1011 uInt k = nrow+i;
1012 scannocol.put(k, newscanno);
1013 const TableRecord& rec = row.get(i);
1014 Double rv,rp,inc;
1015 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
1016 uInt id;
1017 id = out->frequencies().addEntry(rp, rv, inc);
1018 freqidcol.put(k,id);
1019 String name,fname;Double rf;
1020 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
1021 id = out->molecules().addEntry(rf, name, fname);
1022 molidcol.put(k, id);
1023 Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
1024 (*it)->focus().getEntry(fax, ftan, frot, fhand,
1025 fmount,fuser, fxy, fxyp,
1026 rec.asuInt("FOCUS_ID"));
1027 id = out->focus().addEntry(fax, ftan, frot, fhand,
1028 fmount,fuser, fxy, fxyp);
1029 focusidcol.put(k, id);
1030 }
1031 ++freqit;
1032 }
1033 ++newscanno;
1034 ++scanit;
1035 }
1036 ++it;
1037 }
1038 return out;
1039}
1040
1041CountedPtr< Scantable >
1042 STMath::invertPhase( const CountedPtr < Scantable >& in )
1043{
1044 return applyToPol(in, &STPol::invertPhase, Float(0.0));
1045}
1046
1047CountedPtr< Scantable >
1048 STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
1049{
1050 return applyToPol(in, &STPol::rotatePhase, Float(phase));
1051}
1052
1053CountedPtr< Scantable >
1054 STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
1055{
1056 return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
1057}
1058
1059CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
1060 STPol::polOperation fptr,
1061 Float phase )
1062{
1063 CountedPtr< Scantable > out = getScantable(in, false);
1064 Table& tout = out->table();
1065 Block<String> cols(4);
1066 cols[0] = String("SCANNO");
1067 cols[1] = String("BEAMNO");
1068 cols[2] = String("IFNO");
1069 cols[3] = String("CYCLENO");
1070 TableIterator iter(tout, cols);
1071 STPol* stpol = STPol::getPolClass(out->factories_, out->getPolType() );
1072 while (!iter.pastEnd()) {
1073 Table t = iter.table();
1074 ArrayColumn<Float> speccol(t, "SPECTRA");
1075 ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
1076 ScalarColumn<Float> parancol(t, "PARANGLE");
1077 Matrix<Float> pols = speccol.getColumn();
1078 try {
1079 stpol->setSpectra(pols);
1080 Float fang,fhand,parang;
1081 fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
1082 fhand = in->focusTable_.getFeedHand(focidcol(0));
1083 parang = parancol(0);
1084 /// @todo re-enable this
1085 // disable total feed angle to support paralactifying Caswell style
1086 stpol->setPhaseCorrections(parang, -parang, fhand);
1087 (stpol->*fptr)(phase);
1088 speccol.putColumn(stpol->getSpectra());
1089 Matrix<Float> tmp = stpol->getSpectra();
1090 } catch (AipsError& e) {
1091 delete stpol;stpol=0;
1092 throw(e);
1093 }
1094 ++iter;
1095 }
1096 delete stpol;stpol=0;
1097 return out;
1098}
1099
1100CountedPtr< Scantable >
1101 STMath::swapPolarisations( const CountedPtr< Scantable > & in )
1102{
1103 CountedPtr< Scantable > out = getScantable(in, false);
1104 Table& tout = out->table();
1105 Table t0 = tout(tout.col("POLNO") == 0);
1106 Table t1 = tout(tout.col("POLNO") == 1);
1107 if ( t0.nrow() != t1.nrow() )
1108 throw(AipsError("Inconsistent number of polarisations"));
1109 ArrayColumn<Float> speccol0(t0, "SPECTRA");
1110 ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
1111 ArrayColumn<Float> speccol1(t1, "SPECTRA");
1112 ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
1113 Matrix<Float> s0 = speccol0.getColumn();
1114 Matrix<uChar> f0 = flagcol0.getColumn();
1115 speccol0.putColumn(speccol1.getColumn());
1116 flagcol0.putColumn(flagcol1.getColumn());
1117 speccol1.putColumn(s0);
1118 flagcol1.putColumn(f0);
1119 return out;
1120}
1121
1122CountedPtr< Scantable >
1123 STMath::averagePolarisations( const CountedPtr< Scantable > & in,
1124 const std::vector<bool>& mask,
1125 const std::string& weight )
1126{
1127 if (in->getPolType() != "linear" || in->npol() != 2 )
1128 throw(AipsError("averagePolarisations can only be applied to two linear polarisations."));
1129 bool insitu = insitu_;
1130 setInsitu(false);
1131 CountedPtr< Scantable > pols = getScantable(in, false);
1132 setInsitu(insitu);
1133 Table& tout = pols->table();
1134 // give all rows the same POLNO
1135 TableVector<uInt> vec(tout, "POLNO");
1136 vec = 0;
1137 pols->table_.rwKeywordSet().define("nPol", Int(1));
1138 std::vector<CountedPtr<Scantable> > vpols;
1139 vpols.push_back(pols);
1140 CountedPtr< Scantable > out = average(vpols, mask, weight, "NONE");
1141 return out;
1142}
1143
1144CountedPtr< Scantable >
1145 STMath::averageBeams( const CountedPtr< Scantable > & in,
1146 const std::vector<bool>& mask,
1147 const std::string& weight )
1148{
1149 bool insitu = insitu_;
1150 setInsitu(false);
1151 CountedPtr< Scantable > beams = getScantable(in, false);
1152 setInsitu(insitu);
1153 Table& tout = beams->table();
1154 // give all rows the same BEAMNO
1155 TableVector<uInt> vec(tout, "BEAMNO");
1156 vec = 0;
1157 beams->table_.rwKeywordSet().define("nBeam", Int(1));
1158 std::vector<CountedPtr<Scantable> > vbeams;
1159 vbeams.push_back(beams);
1160 CountedPtr< Scantable > out = average(vbeams, mask, weight, "NONE");
1161 return out;
1162}
1163
1164
1165CountedPtr< Scantable >
1166 asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
1167 const std::string & refTime,
1168 const std::string & method)
1169{
1170 // clone as this is not working insitu
1171 bool insitu = insitu_;
1172 setInsitu(false);
1173 CountedPtr< Scantable > out = getScantable(in, false);
1174 setInsitu(insitu);
1175 Table& tout = out->table();
1176 // Get reference Epoch to time of first row or given String
1177 Unit DAY(String("d"));
1178 MEpoch::Ref epochRef(in->getTimeReference());
1179 MEpoch refEpoch;
1180 if (refTime.length()>0) {
1181 Quantum<Double> qt;
1182 if (MVTime::read(qt,refTime)) {
1183 MVEpoch mv(qt);
1184 refEpoch = MEpoch(mv, epochRef);
1185 } else {
1186 throw(AipsError("Invalid format for Epoch string"));
1187 }
1188 } else {
1189 refEpoch = in->timeCol_(0);
1190 }
1191 MPosition refPos = in->getAntennaPosition();
1192
1193 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1194 // test if user frame is different to base frame
1195 if ( in->frequencies().getFrameString(true)
1196 == in->frequencies().getFrameString(false) ) {
1197 throw(AipsError("Can't convert as no output frame has been set"
1198 " (use set_freqframe) or it is aligned already."));
1199 }
1200 MFrequency::Types system = in->frequencies().getFrame();
1201 MVTime mvt(refEpoch.getValue());
1202 String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
1203 ostringstream oss;
1204 oss << "Aligned at reference Epoch " << epochout
1205 << " in frame " << MFrequency::showType(system);
1206 pushLog(String(oss));
1207 // set up the iterator
1208 Block<String> cols(4);
1209 // select by constant direction
1210 cols[0] = String("SRCNAME");
1211 cols[1] = String("BEAMNO");
1212 // select by IF ( no of channels varies over this )
1213 cols[2] = String("IFNO");
1214 // select by restfrequency
1215 cols[3] = String("MOLECULE_ID");
1216 TableIterator iter(tout, cols);
1217 while ( !iter.pastEnd() ) {
1218 Table t = iter.table();
1219 MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1220 TableIterator fiter(t, "FREQ_ID");
1221 // determine nchan from the first row. This should work as
1222 // we are iterating over BEAMNO and IFNO // we should have constant direction
1223
1224 ROArrayColumn<Float> sCol(t, "SPECTRA");
1225 MDirection direction = dirCol(0);
1226 uInt nchan = sCol(0).nelements();
1227 while ( !fiter.pastEnd() ) {
1228 Table ftab = fiter.table();
1229 ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
1230 // get the SpectralCoordinate for the freqid, which we are iterating over
1231 SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1232 FrequencyAligner<Float> fa( sC, nchan, refEpoch,
1233 direction, refPos, system );
1234 // realign the SpectralCoordinate and put into the output Scantable
1235 Vector<String> units(1);
1236 units = String("Hz");
1237 Bool linear=True;
1238 SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1239 sc2.setWorldAxisUnits(units);
1240 uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
1241 sc2.referenceValue()[0],
1242 sc2.increment()[0]);
1243 TableVector<uInt> tvec(ftab, "FREQ_ID");
1244 tvec = id;
1245 // create the "global" abcissa for alignment with same FREQ_ID
1246 Vector<Double> abc(nchan);
1247 Double w;
1248 for (uInt i=0; i<nchan; i++) {
1249 sC.toWorld(w,Double(i));
1250 abc[i] = w;
1251 }
1252 // cache abcissa for same time stamps, so iterate over those
1253 TableIterator timeiter(ftab, "TIME");
1254 while ( !timeiter.pastEnd() ) {
1255 Table tab = timeiter.table();
1256 ArrayColumn<Float> specCol(tab, "SPECTRA");
1257 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1258 MEpoch::ROScalarColumn timeCol(tab, "TIME");
1259 // use align abcissa cache after the first row
1260 bool first = true;
1261 // these rows should be just be POLNO
1262 for (int i=0; i<int(tab.nrow()); ++i) {
1263 // input values
1264 Vector<uChar> flag = flagCol(i);
1265 Vector<Bool> mask(flag.shape());
1266 Vector<Float> specOut, spec;
1267 spec = specCol(i);
1268 Vector<Bool> maskOut;Vector<uChar> flagOut;
1269 convertArray(mask, flag);
1270 // alignment
1271 Bool ok = fa.align(specOut, maskOut, abc, spec,
1272 mask, timeCol(i), !first,
1273 interp, False);
1274 // back into scantable
1275 flagOut.resize(maskOut.nelements());
1276 convertArray(flagOut, maskOut);
1277 flagCol.put(i, flagOut);
1278 specCol.put(i, specOut);
1279 // start abcissa caching
1280 first = false;
1281 }
1282 // next timestamp
1283 ++timeiter;
1284 }
1285 // next FREQ_ID
1286 ++fiter;
1287 }
1288 // next aligner
1289 ++iter;
1290 }
1291 // set this afterwards to ensure we are doing insitu correctly.
1292 out->frequencies().setFrame(system, true);
1293 return out;
1294}
1295
1296CountedPtr<Scantable>
1297 asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
1298 const std::string & newtype )
1299{
1300 if (in->npol() != 2 && in->npol() != 4)
1301 throw(AipsError("Can only convert two or four polarisations."));
1302 if ( in->getPolType() == newtype )
1303 throw(AipsError("No need to convert."));
1304 if ( ! in->selector_.empty() )
1305 throw(AipsError("Can only convert whole scantable. Unset the selection."));
1306 bool insitu = insitu_;
1307 setInsitu(false);
1308 CountedPtr< Scantable > out = getScantable(in, true);
1309 setInsitu(insitu);
1310 Table& tout = out->table();
1311 tout.rwKeywordSet().define("POLTYPE", String(newtype));
1312
1313 Block<String> cols(4);
1314 cols[0] = "SCANNO";
1315 cols[1] = "CYCLENO";
1316 cols[2] = "BEAMNO";
1317 cols[3] = "IFNO";
1318 TableIterator it(in->originalTable_, cols);
1319 String basetype = in->getPolType();
1320 STPol* stpol = STPol::getPolClass(in->factories_, basetype);
1321 try {
1322 while ( !it.pastEnd() ) {
1323 Table tab = it.table();
1324 uInt row = tab.rowNumbers()[0];
1325 stpol->setSpectra(in->getPolMatrix(row));
1326 Float fang,fhand,parang;
1327 fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
1328 fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
1329 parang = in->paraCol_(row);
1330 /// @todo re-enable this
1331 // disable total feed angle to support paralactifying Caswell style
1332 stpol->setPhaseCorrections(parang, -parang, fhand);
1333 Int npolout = 0;
1334 for (uInt i=0; i<tab.nrow(); ++i) {
1335 Vector<Float> outvec = stpol->getSpectrum(i, newtype);
1336 if ( outvec.nelements() > 0 ) {
1337 tout.addRow();
1338 TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
1339 ArrayColumn<Float> sCol(tout,"SPECTRA");
1340 ScalarColumn<uInt> pCol(tout,"POLNO");
1341 sCol.put(tout.nrow()-1 ,outvec);
1342 pCol.put(tout.nrow()-1 ,uInt(npolout));
1343 npolout++;
1344 }
1345 }
1346 tout.rwKeywordSet().define("nPol", npolout);
1347 ++it;
1348 }
1349 } catch (AipsError& e) {
1350 delete stpol;
1351 throw(e);
1352 }
1353 delete stpol;
1354 return out;
1355}
1356
1357CountedPtr< Scantable >
1358 asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
1359 const std::string & scantype )
1360{
1361 bool insitu = insitu_;
1362 setInsitu(false);
1363 CountedPtr< Scantable > out = getScantable(in, true);
1364 setInsitu(insitu);
1365 Table& tout = out->table();
1366 std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
1367 if (scantype == "on") {
1368 taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
1369 }
1370 Table tab = tableCommand(taql, in->table());
1371 TableCopy::copyRows(tout, tab);
1372 if (scantype == "on") {
1373 // re-index SCANNO to 0
1374 TableVector<uInt> vec(tout, "SCANNO");
1375 vec = 0;
1376 }
1377 return out;
1378}
Note: See TracBrowser for help on using the repository browser.