source: trunk/src/STMath.cpp@ 1141

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

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