source: branches/Release2.1/src/STMath.cpp@ 2793

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

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