source: trunk/src/STMath.cpp@ 1217

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

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