source: trunk/src/STMath.cpp@ 1069

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

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