source: trunk/src/STMath.cpp@ 1079

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

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