source: trunk/src/STMath.cpp@ 1460

Last change on this file since 1460 was 1439, checked in by Malte Marquarding, 16 years ago

added Header diff function. Don't use obstype for conformance test
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 64.9 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#include "STSelector.h"
51
52using namespace casa;
53
54using namespace asap;
55
56STMath::STMath(bool insitu) :
57 insitu_(insitu)
58{
59}
60
61
62STMath::~STMath()
63{
64}
65
66CountedPtr<Scantable>
67STMath::average( const std::vector<CountedPtr<Scantable> >& in,
68 const std::vector<bool>& mask,
69 const std::string& weight,
70 const std::string& avmode)
71{
72 if ( avmode == "SCAN" && in.size() != 1 )
73 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
74 "Use merge first."));
75 WeightType wtype = stringToWeight(weight);
76
77 // output
78 // clone as this is non insitu
79 bool insitu = insitu_;
80 setInsitu(false);
81 CountedPtr< Scantable > out = getScantable(in[0], true);
82 setInsitu(insitu);
83 std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
84 ++stit;
85 while ( stit != in.end() ) {
86 out->appendToHistoryTable((*stit)->history());
87 ++stit;
88 }
89
90 Table& tout = out->table();
91
92 /// @todo check if all scantables are conformant
93
94 ArrayColumn<Float> specColOut(tout,"SPECTRA");
95 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
96 ArrayColumn<Float> tsysColOut(tout,"TSYS");
97 ScalarColumn<Double> mjdColOut(tout,"TIME");
98 ScalarColumn<Double> intColOut(tout,"INTERVAL");
99 ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
100 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
101
102 // set up the output table rows. These are based on the structure of the
103 // FIRST scantable in the vector
104 const Table& baset = in[0]->table();
105
106 Block<String> cols(3);
107 cols[0] = String("BEAMNO");
108 cols[1] = String("IFNO");
109 cols[2] = String("POLNO");
110 if ( avmode == "SOURCE" ) {
111 cols.resize(4);
112 cols[3] = String("SRCNAME");
113 }
114 if ( avmode == "SCAN" && in.size() == 1) {
115 cols.resize(4);
116 cols[3] = String("SCANNO");
117 }
118 uInt outrowCount = 0;
119 TableIterator iter(baset, cols);
120 while (!iter.pastEnd()) {
121 Table subt = iter.table();
122 // copy the first row of this selection into the new table
123 tout.addRow();
124 TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
125 // re-index to 0
126 if ( avmode != "SCAN" && avmode != "SOURCE" ) {
127 scanColOut.put(outrowCount, uInt(0));
128 }
129 ++outrowCount;
130 ++iter;
131 }
132 RowAccumulator acc(wtype);
133 Vector<Bool> cmask(mask);
134 acc.setUserMask(cmask);
135 ROTableRow row(tout);
136 ROArrayColumn<Float> specCol, tsysCol;
137 ROArrayColumn<uChar> flagCol;
138 ROScalarColumn<Double> mjdCol, intCol;
139 ROScalarColumn<Int> scanIDCol;
140
141 Vector<uInt> rowstodelete;
142
143 for (uInt i=0; i < tout.nrow(); ++i) {
144 for ( int j=0; j < int(in.size()); ++j ) {
145 const Table& tin = in[j]->table();
146 const TableRecord& rec = row.get(i);
147 ROScalarColumn<Double> tmp(tin, "TIME");
148 Double td;tmp.get(0,td);
149 Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
150 && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
151 && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
152 Table subt;
153 if ( avmode == "SOURCE") {
154 subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
155 } else if (avmode == "SCAN") {
156 subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
157 } else {
158 subt = basesubt;
159 }
160 specCol.attach(subt,"SPECTRA");
161 flagCol.attach(subt,"FLAGTRA");
162 tsysCol.attach(subt,"TSYS");
163 intCol.attach(subt,"INTERVAL");
164 mjdCol.attach(subt,"TIME");
165 Vector<Float> spec,tsys;
166 Vector<uChar> flag;
167 Double inter,time;
168 for (uInt k = 0; k < subt.nrow(); ++k ) {
169 flagCol.get(k, flag);
170 Vector<Bool> bflag(flag.shape());
171 convertArray(bflag, flag);
172 /*
173 if ( allEQ(bflag, True) ) {
174 continue;//don't accumulate
175 }
176 */
177 specCol.get(k, spec);
178 tsysCol.get(k, tsys);
179 intCol.get(k, inter);
180 mjdCol.get(k, time);
181 // spectrum has to be added last to enable weighting by the other values
182 acc.add(spec, !bflag, tsys, inter, time);
183 }
184 }
185 const Vector<Bool>& msk = acc.getMask();
186 if ( allEQ(msk, False) ) {
187 uint n = rowstodelete.nelements();
188 rowstodelete.resize(n+1, True);
189 rowstodelete[n] = i;
190 continue;
191 }
192 //write out
193 Vector<uChar> flg(msk.shape());
194 convertArray(flg, !msk);
195 flagColOut.put(i, flg);
196 specColOut.put(i, acc.getSpectrum());
197 tsysColOut.put(i, acc.getTsys());
198 intColOut.put(i, acc.getInterval());
199 mjdColOut.put(i, acc.getTime());
200 // we should only have one cycle now -> reset it to be 0
201 // frequency switched data has different CYCLENO for different IFNO
202 // which requires resetting this value
203 cycColOut.put(i, uInt(0));
204 acc.reset();
205 }
206 if (rowstodelete.nelements() > 0) {
207 tout.removeRow(rowstodelete);
208 if (tout.nrow() == 0) {
209 throw(AipsError("Can't average fully flagged data."));
210 }
211 }
212 return out;
213}
214
215CountedPtr< Scantable >
216 STMath::averageChannel( const CountedPtr < Scantable > & in,
217 const std::string & mode,
218 const std::string& avmode )
219{
220 // clone as this is non insitu
221 bool insitu = insitu_;
222 setInsitu(false);
223 CountedPtr< Scantable > out = getScantable(in, true);
224 setInsitu(insitu);
225 Table& tout = out->table();
226 ArrayColumn<Float> specColOut(tout,"SPECTRA");
227 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
228 ArrayColumn<Float> tsysColOut(tout,"TSYS");
229 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
230 ScalarColumn<Double> intColOut(tout, "INTERVAL");
231 Table tmp = in->table().sort("BEAMNO");
232 Block<String> cols(3);
233 cols[0] = String("BEAMNO");
234 cols[1] = String("IFNO");
235 cols[2] = String("POLNO");
236 if ( avmode == "SCAN") {
237 cols.resize(4);
238 cols[3] = String("SCANNO");
239 }
240 uInt outrowCount = 0;
241 uChar userflag = 1 << 7;
242 TableIterator iter(tmp, cols);
243 while (!iter.pastEnd()) {
244 Table subt = iter.table();
245 ROArrayColumn<Float> specCol, tsysCol;
246 ROArrayColumn<uChar> flagCol;
247 ROScalarColumn<Double> intCol(subt, "INTERVAL");
248 specCol.attach(subt,"SPECTRA");
249 flagCol.attach(subt,"FLAGTRA");
250 tsysCol.attach(subt,"TSYS");
251 tout.addRow();
252 TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
253 if ( avmode != "SCAN") {
254 scanColOut.put(outrowCount, uInt(0));
255 }
256 Vector<Float> tmp;
257 specCol.get(0, tmp);
258 uInt nchan = tmp.nelements();
259 // have to do channel by channel here as MaskedArrMath
260 // doesn't have partialMedians
261 Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
262 Vector<Float> outspec(nchan);
263 Vector<uChar> outflag(nchan,0);
264 Vector<Float> outtsys(1);/// @fixme when tsys is channel based
265 for (uInt i=0; i<nchan; ++i) {
266 Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
267 MaskedArray<Float> ma = maskedArray(specs,flags);
268 outspec[i] = median(ma);
269 if ( allEQ(ma.getMask(), False) )
270 outflag[i] = userflag;// flag data
271 }
272 outtsys[0] = median(tsysCol.getColumn());
273 specColOut.put(outrowCount, outspec);
274 flagColOut.put(outrowCount, outflag);
275 tsysColOut.put(outrowCount, outtsys);
276 Double intsum = sum(intCol.getColumn());
277 intColOut.put(outrowCount, intsum);
278 ++outrowCount;
279 ++iter;
280 }
281 return out;
282}
283
284CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
285 bool droprows)
286{
287 if (insitu_) return in;
288 else {
289 // clone
290 Scantable* tabp = new Scantable(*in, Bool(droprows));
291 return CountedPtr<Scantable>(tabp);
292 }
293}
294
295CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
296 float val,
297 const std::string& mode,
298 bool tsys )
299{
300 CountedPtr< Scantable > out = getScantable(in, false);
301 Table& tab = out->table();
302 ArrayColumn<Float> specCol(tab,"SPECTRA");
303 ArrayColumn<Float> tsysCol(tab,"TSYS");
304 for (uInt i=0; i<tab.nrow(); ++i) {
305 Vector<Float> spec;
306 Vector<Float> ts;
307 specCol.get(i, spec);
308 tsysCol.get(i, ts);
309 if (mode == "MUL" || mode == "DIV") {
310 if (mode == "DIV") val = 1.0/val;
311 spec *= val;
312 specCol.put(i, spec);
313 if ( tsys ) {
314 ts *= val;
315 tsysCol.put(i, ts);
316 }
317 } else if ( mode == "ADD" || mode == "SUB") {
318 if (mode == "SUB") val *= -1.0;
319 spec += val;
320 specCol.put(i, spec);
321 if ( tsys ) {
322 ts += val;
323 tsysCol.put(i, ts);
324 }
325 }
326 }
327 return out;
328}
329
330CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
331 const CountedPtr<Scantable>& right,
332 const std::string& mode)
333{
334 bool insitu = insitu_;
335 if ( ! left->conformant(*right) ) {
336 throw(AipsError("'left' and 'right' scantables are not conformant."));
337 }
338 setInsitu(false);
339 CountedPtr< Scantable > out = getScantable(left, false);
340 setInsitu(insitu);
341 Table& tout = out->table();
342 Block<String> coln(5);
343 coln[0] = "SCANNO"; coln[1] = "CYCLENO"; coln[2] = "BEAMNO";
344 coln[3] = "IFNO"; coln[4] = "POLNO";
345 Table tmpl = tout.sort(coln);
346 Table tmpr = right->table().sort(coln);
347 ArrayColumn<Float> lspecCol(tmpl,"SPECTRA");
348 ROArrayColumn<Float> rspecCol(tmpr,"SPECTRA");
349 ArrayColumn<uChar> lflagCol(tmpl,"FLAGTRA");
350 ROArrayColumn<uChar> rflagCol(tmpr,"FLAGTRA");
351
352 for (uInt i=0; i<tout.nrow(); ++i) {
353 Vector<Float> lspecvec, rspecvec;
354 Vector<uChar> lflagvec, rflagvec;
355 lspecvec = lspecCol(i); rspecvec = rspecCol(i);
356 lflagvec = lflagCol(i); rflagvec = rflagCol(i);
357 MaskedArray<Float> mleft = maskedArray(lspecvec, lflagvec);
358 MaskedArray<Float> mright = maskedArray(rspecvec, rflagvec);
359 if (mode == "ADD") {
360 mleft += mright;
361 } else if ( mode == "SUB") {
362 mleft -= mright;
363 } else if ( mode == "MUL") {
364 mleft *= mright;
365 } else if ( mode == "DIV") {
366 mleft /= mright;
367 } else {
368 throw(AipsError("Illegal binary operator"));
369 }
370 lspecCol.put(i, mleft.getArray());
371 }
372 return out;
373}
374
375
376
377MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
378 const Vector<uChar>& f)
379{
380 Vector<Bool> mask;
381 mask.resize(f.shape());
382 convertArray(mask, f);
383 return MaskedArray<Float>(s,!mask);
384}
385
386Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
387{
388 const Vector<Bool>& m = ma.getMask();
389 Vector<uChar> flags(m.shape());
390 convertArray(flags, !m);
391 return flags;
392}
393
394CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
395 const std::string & mode,
396 bool preserve )
397{
398 /// @todo make other modes available
399 /// modes should be "nearest", "pair"
400 // make this operation non insitu
401 const Table& tin = in->table();
402 Table ons = tin(tin.col("SRCTYPE") == Int(0));
403 Table offs = tin(tin.col("SRCTYPE") == Int(1));
404 if ( offs.nrow() == 0 )
405 throw(AipsError("No 'off' scans present."));
406 // put all "on" scans into output table
407
408 bool insitu = insitu_;
409 setInsitu(false);
410 CountedPtr< Scantable > out = getScantable(in, true);
411 setInsitu(insitu);
412 Table& tout = out->table();
413
414 TableCopy::copyRows(tout, ons);
415 TableRow row(tout);
416 ROScalarColumn<Double> offtimeCol(offs, "TIME");
417 ArrayColumn<Float> outspecCol(tout, "SPECTRA");
418 ROArrayColumn<Float> outtsysCol(tout, "TSYS");
419 ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
420 for (uInt i=0; i < tout.nrow(); ++i) {
421 const TableRecord& rec = row.get(i);
422 Double ontime = rec.asDouble("TIME");
423 Table presel = offs(offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
424 && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
425 && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
426 ROScalarColumn<Double> offtimeCol(presel, "TIME");
427
428 Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
429 // Timestamp may vary within a cycle ???!!!
430 // increase this by 0.01 sec in case of rounding errors...
431 // There might be a better way to do this.
432 // fix to this fix. TIME is MJD, so 1.0d not 1.0s
433 mindeltat += 0.01/24./60./60.;
434 Table sel = presel( abs(presel.col("TIME")-ontime) <= mindeltat);
435
436 if ( sel.nrow() < 1 ) {
437 throw(AipsError("No closest in time found... This could be a rounding "
438 "issue. Try quotient instead."));
439 }
440 TableRow offrow(sel);
441 const TableRecord& offrec = offrow.get(0);//should only be one row
442 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
443 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
444 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
445 /// @fixme this assumes tsys is a scalar not vector
446 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
447 Vector<Float> specon, tsyson;
448 outtsysCol.get(i, tsyson);
449 outspecCol.get(i, specon);
450 Vector<uChar> flagon;
451 outflagCol.get(i, flagon);
452 MaskedArray<Float> mon = maskedArray(specon, flagon);
453 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
454 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
455 if (preserve) {
456 quot -= tsysoffscalar;
457 } else {
458 quot -= tsyson[0];
459 }
460 outspecCol.put(i, quot.getArray());
461 outflagCol.put(i, flagsFromMA(quot));
462 }
463 // renumber scanno
464 TableIterator it(tout, "SCANNO");
465 uInt i = 0;
466 while ( !it.pastEnd() ) {
467 Table t = it.table();
468 TableVector<uInt> vec(t, "SCANNO");
469 vec = i;
470 ++i;
471 ++it;
472 }
473 return out;
474}
475
476
477CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
478 const CountedPtr< Scantable > & off,
479 bool preserve )
480{
481 bool insitu = insitu_;
482 if ( ! on->conformant(*off) ) {
483 throw(AipsError("'on' and 'off' scantables are not conformant."));
484 }
485 setInsitu(false);
486 CountedPtr< Scantable > out = getScantable(on, false);
487 setInsitu(insitu);
488 Table& tout = out->table();
489 const Table& toff = off->table();
490 TableIterator sit(tout, "SCANNO");
491 TableIterator s2it(toff, "SCANNO");
492 while ( !sit.pastEnd() ) {
493 Table ton = sit.table();
494 TableRow row(ton);
495 Table t = s2it.table();
496 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
497 ROArrayColumn<Float> outtsysCol(ton, "TSYS");
498 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
499 for (uInt i=0; i < ton.nrow(); ++i) {
500 const TableRecord& rec = row.get(i);
501 Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
502 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
503 && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
504 if ( offsel.nrow() == 0 )
505 throw AipsError("STMath::quotient: no matching off");
506 TableRow offrow(offsel);
507 const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
508 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
509 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
510 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
511 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
512 Vector<Float> specon, tsyson;
513 outtsysCol.get(i, tsyson);
514 outspecCol.get(i, specon);
515 Vector<uChar> flagon;
516 outflagCol.get(i, flagon);
517 MaskedArray<Float> mon = maskedArray(specon, flagon);
518 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
519 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
520 if (preserve) {
521 quot -= tsysoffscalar;
522 } else {
523 quot -= tsyson[0];
524 }
525 outspecCol.put(i, quot.getArray());
526 outflagCol.put(i, flagsFromMA(quot));
527 }
528 ++sit;
529 ++s2it;
530 // take the first off for each on scan which doesn't have a
531 // matching off scan
532 // non <= noff: matching pairs, non > noff matching pairs then first off
533 if ( s2it.pastEnd() ) s2it.reset();
534 }
535 return out;
536}
537
538// dototalpower (migration of GBTIDL procedure dototalpower.pro)
539// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
540// do it for each cycles in a specific scan.
541CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
542 const CountedPtr< Scantable >& caloff, Float tcal )
543{
544if ( ! calon->conformant(*caloff) ) {
545 throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
546 }
547 setInsitu(false);
548 CountedPtr< Scantable > out = getScantable(caloff, false);
549 Table& tout = out->table();
550 const Table& tcon = calon->table();
551 Vector<Float> tcalout;
552 Vector<Float> tcalout2; //debug
553
554 if ( tout.nrow() != tcon.nrow() ) {
555 throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
556 }
557 // iteration by scanno or cycle no.
558 TableIterator sit(tout, "SCANNO");
559 TableIterator s2it(tcon, "SCANNO");
560 while ( !sit.pastEnd() ) {
561 Table toff = sit.table();
562 TableRow row(toff);
563 Table t = s2it.table();
564 ScalarColumn<Double> outintCol(toff, "INTERVAL");
565 ArrayColumn<Float> outspecCol(toff, "SPECTRA");
566 ArrayColumn<Float> outtsysCol(toff, "TSYS");
567 ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
568 ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
569 ROScalarColumn<uInt> outpolCol(toff, "POLNO");
570 ROScalarColumn<Double> onintCol(t, "INTERVAL");
571 ROArrayColumn<Float> onspecCol(t, "SPECTRA");
572 ROArrayColumn<Float> ontsysCol(t, "TSYS");
573 ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
574 //ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
575
576 for (uInt i=0; i < toff.nrow(); ++i) {
577 //skip these checks -> assumes the data order are the same between the cal on off pairs
578 //
579 Vector<Float> specCalon, specCaloff;
580 // to store scalar (mean) tsys
581 Vector<Float> tsysout(1);
582 uInt tcalId, polno;
583 Double offint, onint;
584 outpolCol.get(i, polno);
585 outspecCol.get(i, specCaloff);
586 onspecCol.get(i, specCalon);
587 Vector<uChar> flagCaloff, flagCalon;
588 outflagCol.get(i, flagCaloff);
589 onflagCol.get(i, flagCalon);
590 outtcalIdCol.get(i, tcalId);
591 outintCol.get(i, offint);
592 onintCol.get(i, onint);
593 // caluculate mean Tsys
594 uInt nchan = specCaloff.nelements();
595 // percentage of edge cut off
596 uInt pc = 10;
597 uInt bchan = nchan/pc;
598 uInt echan = nchan-bchan;
599
600 Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
601 Vector<Float> testsubsp = specCaloff(chansl);
602 MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
603 MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
604 MaskedArray<Float> spdiff = spon-spoff;
605 uInt noff = spoff.nelementsValid();
606 //uInt non = spon.nelementsValid();
607 uInt ndiff = spdiff.nelementsValid();
608 Float meantsys;
609
610/**
611 Double subspec, subdiff;
612 uInt usednchan;
613 subspec = 0;
614 subdiff = 0;
615 usednchan = 0;
616 for(uInt k=(bchan-1); k<echan; k++) {
617 subspec += specCaloff[k];
618 subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
619 ++usednchan;
620 }
621**/
622 // get tcal if input tcal <= 0
623 String tcalt;
624 Float tcalUsed;
625 tcalUsed = tcal;
626 if ( tcal <= 0.0 ) {
627 caloff->tcal().getEntry(tcalt, tcalout, tcalId);
628 if (polno<=3) {
629 tcalUsed = tcalout[polno];
630 }
631 else {
632 tcalUsed = tcalout[0];
633 }
634 }
635
636 Float meanoff;
637 Float meandiff;
638 if (noff && ndiff) {
639 //Debug
640 //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
641 meanoff = sum(spoff)/noff;
642 meandiff = sum(spdiff)/ndiff;
643 meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
644 }
645 else {
646 meantsys=1;
647 }
648
649 tsysout[0] = Float(meantsys);
650 MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
651 MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
652 MaskedArray<Float> sig = Float(0.5) * (mcaloff + mcalon);
653 //uInt ncaloff = mcaloff.nelementsValid();
654 //uInt ncalon = mcalon.nelementsValid();
655
656 outintCol.put(i, offint+onint);
657 outspecCol.put(i, sig.getArray());
658 outflagCol.put(i, flagsFromMA(sig));
659 outtsysCol.put(i, tsysout);
660 }
661 ++sit;
662 ++s2it;
663 }
664 return out;
665}
666
667//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
668// observatiions.
669// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
670// no smoothing).
671// output: resultant scantable [= (sig-ref/ref)*tsys]
672CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
673 const CountedPtr < Scantable >& ref,
674 int smoothref,
675 casa::Float tsysv,
676 casa::Float tau )
677{
678if ( ! ref->conformant(*sig) ) {
679 throw(AipsError("'sig' and 'ref' scantables are not conformant."));
680 }
681 setInsitu(false);
682 CountedPtr< Scantable > out = getScantable(sig, false);
683 CountedPtr< Scantable > smref;
684 if ( smoothref > 1 ) {
685 float fsmoothref = static_cast<float>(smoothref);
686 std::string inkernel = "boxcar";
687 smref = smooth(ref, inkernel, fsmoothref );
688 ostringstream oss;
689 oss<<"Applied smoothing of "<<fsmoothref<<" on the reference."<<endl;
690 pushLog(String(oss));
691 }
692 else {
693 smref = ref;
694 }
695 Table& tout = out->table();
696 const Table& tref = smref->table();
697 if ( tout.nrow() != tref.nrow() ) {
698 throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
699 }
700 // iteration by scanno? or cycle no.
701 TableIterator sit(tout, "SCANNO");
702 TableIterator s2it(tref, "SCANNO");
703 while ( !sit.pastEnd() ) {
704 Table ton = sit.table();
705 Table t = s2it.table();
706 ScalarColumn<Double> outintCol(ton, "INTERVAL");
707 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
708 ArrayColumn<Float> outtsysCol(ton, "TSYS");
709 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
710 ArrayColumn<Float> refspecCol(t, "SPECTRA");
711 ROScalarColumn<Double> refintCol(t, "INTERVAL");
712 ROArrayColumn<Float> reftsysCol(t, "TSYS");
713 ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
714 ROScalarColumn<Float> refelevCol(t, "ELEVATION");
715 for (uInt i=0; i < ton.nrow(); ++i) {
716
717 Double onint, refint;
718 Vector<Float> specon, specref;
719 // to store scalar (mean) tsys
720 Vector<Float> tsysref;
721 outintCol.get(i, onint);
722 refintCol.get(i, refint);
723 outspecCol.get(i, specon);
724 refspecCol.get(i, specref);
725 Vector<uChar> flagref, flagon;
726 outflagCol.get(i, flagon);
727 refflagCol.get(i, flagref);
728 reftsysCol.get(i, tsysref);
729
730 Float tsysrefscalar;
731 if ( tsysv > 0.0 ) {
732 ostringstream oss;
733 Float elev;
734 refelevCol.get(i, elev);
735 oss << "user specified Tsys = " << tsysv;
736 // do recalc elevation if EL = 0
737 if ( elev == 0 ) {
738 throw(AipsError("EL=0, elevation data is missing."));
739 } else {
740 if ( tau <= 0.0 ) {
741 throw(AipsError("Valid tau is not supplied."));
742 } else {
743 tsysrefscalar = tsysv * exp(tau/elev);
744 }
745 }
746 oss << ", corrected (for El) tsys= "<<tsysrefscalar;
747 pushLog(String(oss));
748 }
749 else {
750 tsysrefscalar = tsysref[0];
751 }
752 //get quotient spectrum
753 MaskedArray<Float> mref = maskedArray(specref, flagref);
754 MaskedArray<Float> mon = maskedArray(specon, flagon);
755 MaskedArray<Float> specres = tsysrefscalar*((mon - mref)/mref);
756 Double resint = onint*refint*smoothref/(onint+refint*smoothref);
757
758 //Debug
759 //cerr<<"Tsys used="<<tsysrefscalar<<endl;
760 // fill the result, replay signal tsys by reference tsys
761 outintCol.put(i, resint);
762 outspecCol.put(i, specres.getArray());
763 outflagCol.put(i, flagsFromMA(specres));
764 outtsysCol.put(i, tsysref);
765 }
766 ++sit;
767 ++s2it;
768 }
769 return out;
770}
771
772CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
773 const std::vector<int>& scans,
774 int smoothref,
775 casa::Float tsysv,
776 casa::Float tau,
777 casa::Float tcal )
778
779{
780 setInsitu(false);
781 STSelector sel;
782 std::vector<int> scan1, scan2, beams;
783 std::vector< vector<int> > scanpair;
784 std::vector<string> calstate;
785 String msg;
786
787 CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
788 CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
789
790 std::vector< CountedPtr< Scantable > > sctables;
791 sctables.push_back(s1b1on);
792 sctables.push_back(s1b1off);
793 sctables.push_back(s1b2on);
794 sctables.push_back(s1b2off);
795 sctables.push_back(s2b1on);
796 sctables.push_back(s2b1off);
797 sctables.push_back(s2b2on);
798 sctables.push_back(s2b2off);
799
800 //check scanlist
801 int n=s->checkScanInfo(scans);
802 if (n==1) {
803 throw(AipsError("Incorrect scan pairs. "));
804 }
805
806 // Assume scans contain only a pair of consecutive scan numbers.
807 // It is assumed that first beam, b1, is on target.
808 // There is no check if the first beam is on or not.
809 if ( scans.size()==1 ) {
810 scan1.push_back(scans[0]);
811 scan2.push_back(scans[0]+1);
812 } else if ( scans.size()==2 ) {
813 scan1.push_back(scans[0]);
814 scan2.push_back(scans[1]);
815 } else {
816 if ( scans.size()%2 == 0 ) {
817 for (uInt i=0; i<scans.size(); i++) {
818 if (i%2 == 0) {
819 scan1.push_back(scans[i]);
820 }
821 else {
822 scan2.push_back(scans[i]);
823 }
824 }
825 } else {
826 throw(AipsError("Odd numbers of scans, cannot form pairs."));
827 }
828 }
829 scanpair.push_back(scan1);
830 scanpair.push_back(scan2);
831 calstate.push_back("*calon");
832 calstate.push_back("*[^calon]");
833 CountedPtr< Scantable > ws = getScantable(s, false);
834 uInt l=0;
835 while ( l < sctables.size() ) {
836 for (uInt i=0; i < 2; i++) {
837 for (uInt j=0; j < 2; j++) {
838 for (uInt k=0; k < 2; k++) {
839 sel.reset();
840 sel.setScans(scanpair[i]);
841 sel.setName(calstate[k]);
842 beams.clear();
843 beams.push_back(j);
844 sel.setBeams(beams);
845 ws->setSelection(sel);
846 sctables[l]= getScantable(ws, false);
847 l++;
848 }
849 }
850 }
851 }
852
853 // replace here by splitData or getData functionality
854 CountedPtr< Scantable > sig1;
855 CountedPtr< Scantable > ref1;
856 CountedPtr< Scantable > sig2;
857 CountedPtr< Scantable > ref2;
858 CountedPtr< Scantable > calb1;
859 CountedPtr< Scantable > calb2;
860
861 msg=String("Processing dototalpower for subset of the data");
862 ostringstream oss1;
863 oss1 << msg << endl;
864 pushLog(String(oss1));
865 // Debug for IRC CS data
866 //float tcal1=7.0;
867 //float tcal2=4.0;
868 sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
869 ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
870 ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
871 sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
872
873 // correction of user-specified tsys for elevation here
874
875 // dosigref calibration
876 msg=String("Processing dosigref for subset of the data");
877 ostringstream oss2;
878 oss2 << msg << endl;
879 pushLog(String(oss2));
880 calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
881 calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
882
883 // iteration by scanno or cycle no.
884 Table& tcalb1 = calb1->table();
885 Table& tcalb2 = calb2->table();
886 TableIterator sit(tcalb1, "SCANNO");
887 TableIterator s2it(tcalb2, "SCANNO");
888 while ( !sit.pastEnd() ) {
889 Table t1 = sit.table();
890 Table t2= s2it.table();
891 ArrayColumn<Float> outspecCol(t1, "SPECTRA");
892 ArrayColumn<Float> outtsysCol(t1, "TSYS");
893 ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
894 ScalarColumn<Double> outintCol(t1, "INTERVAL");
895 ArrayColumn<Float> t2specCol(t2, "SPECTRA");
896 ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
897 ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
898 ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
899 for (uInt i=0; i < t1.nrow(); ++i) {
900 Vector<Float> spec1, spec2;
901 // to store scalar (mean) tsys
902 Vector<Float> tsys1, tsys2;
903 Vector<uChar> flag1, flag2;
904 Double tint1, tint2;
905 outspecCol.get(i, spec1);
906 t2specCol.get(i, spec2);
907 outflagCol.get(i, flag1);
908 t2flagCol.get(i, flag2);
909 outtsysCol.get(i, tsys1);
910 t2tsysCol.get(i, tsys2);
911 outintCol.get(i, tint1);
912 t2intCol.get(i, tint2);
913 // average
914 // assume scalar tsys for weights
915 Float wt1, wt2, tsyssq1, tsyssq2;
916 tsyssq1 = tsys1[0]*tsys1[0];
917 tsyssq2 = tsys2[0]*tsys2[0];
918 wt1 = Float(tint1)/tsyssq1;
919 wt2 = Float(tint2)/tsyssq2;
920 Float invsumwt=1/(wt1+wt2);
921 MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
922 MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
923 MaskedArray<Float> avspec = invsumwt * (wt1*mspec1 + wt2*mspec2);
924 //Array<Float> avtsys = Float(0.5) * (tsys1 + tsys2);
925 // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
926 tsys1[0] = sqrt(tsyssq1 + tsyssq2);
927 Array<Float> avtsys = tsys1;
928
929 outspecCol.put(i, avspec.getArray());
930 outflagCol.put(i, flagsFromMA(avspec));
931 outtsysCol.put(i, avtsys);
932 }
933 ++sit;
934 ++s2it;
935 }
936 return calb1;
937}
938
939//GBTIDL version of frequency switched data calibration
940CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
941 const std::vector<int>& scans,
942 int smoothref,
943 casa::Float tsysv,
944 casa::Float tau,
945 casa::Float tcal )
946{
947
948
949 STSelector sel;
950 CountedPtr< Scantable > ws = getScantable(s, false);
951 CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
952 CountedPtr< Scantable > calsig, calref, out;
953
954 //split the data
955 sel.setName("*_fs");
956 ws->setSelection(sel);
957 sig = getScantable(ws,false);
958 sel.reset();
959 sel.setName("*_fs_calon");
960 ws->setSelection(sel);
961 sigwcal = getScantable(ws,false);
962 sel.reset();
963 sel.setName("*_fsr");
964 ws->setSelection(sel);
965 ref = getScantable(ws,false);
966 sel.reset();
967 sel.setName("*_fsr_calon");
968 ws->setSelection(sel);
969 refwcal = getScantable(ws,false);
970
971 calsig = dototalpower(sigwcal, sig, tcal=tcal);
972 calref = dototalpower(refwcal, ref, tcal=tcal);
973
974 out=dosigref(calsig,calref,smoothref,tsysv,tau);
975
976 return out;
977}
978
979
980CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
981{
982 // make copy or reference
983 CountedPtr< Scantable > out = getScantable(in, false);
984 Table& tout = out->table();
985 Block<String> cols(4);
986 cols[0] = String("SCANNO");
987 cols[1] = String("CYCLENO");
988 cols[2] = String("BEAMNO");
989 cols[3] = String("POLNO");
990 TableIterator iter(tout, cols);
991 while (!iter.pastEnd()) {
992 Table subt = iter.table();
993 // this should leave us with two rows for the two IFs....if not ignore
994 if (subt.nrow() != 2 ) {
995 continue;
996 }
997 ArrayColumn<Float> specCol(subt, "SPECTRA");
998 ArrayColumn<Float> tsysCol(subt, "TSYS");
999 ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
1000 Vector<Float> onspec,offspec, ontsys, offtsys;
1001 Vector<uChar> onflag, offflag;
1002 tsysCol.get(0, ontsys); tsysCol.get(1, offtsys);
1003 specCol.get(0, onspec); specCol.get(1, offspec);
1004 flagCol.get(0, onflag); flagCol.get(1, offflag);
1005 MaskedArray<Float> on = maskedArray(onspec, onflag);
1006 MaskedArray<Float> off = maskedArray(offspec, offflag);
1007 MaskedArray<Float> oncopy = on.copy();
1008
1009 on /= off; on -= 1.0f;
1010 on *= ontsys[0];
1011 off /= oncopy; off -= 1.0f;
1012 off *= offtsys[0];
1013 specCol.put(0, on.getArray());
1014 const Vector<Bool>& m0 = on.getMask();
1015 Vector<uChar> flags0(m0.shape());
1016 convertArray(flags0, !m0);
1017 flagCol.put(0, flags0);
1018
1019 specCol.put(1, off.getArray());
1020 const Vector<Bool>& m1 = off.getMask();
1021 Vector<uChar> flags1(m1.shape());
1022 convertArray(flags1, !m1);
1023 flagCol.put(1, flags1);
1024 ++iter;
1025 }
1026
1027 return out;
1028}
1029
1030std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
1031 const std::vector< bool > & mask,
1032 const std::string& which )
1033{
1034
1035 Vector<Bool> m(mask);
1036 const Table& tab = in->table();
1037 ROArrayColumn<Float> specCol(tab, "SPECTRA");
1038 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1039 std::vector<float> out;
1040 for (uInt i=0; i < tab.nrow(); ++i ) {
1041 Vector<Float> spec; specCol.get(i, spec);
1042 Vector<uChar> flag; flagCol.get(i, flag);
1043 MaskedArray<Float> ma = maskedArray(spec, flag);
1044 float outstat = 0.0;
1045 if ( spec.nelements() == m.nelements() ) {
1046 outstat = mathutil::statistics(which, ma(m));
1047 } else {
1048 outstat = mathutil::statistics(which, ma);
1049 }
1050 out.push_back(outstat);
1051 }
1052 return out;
1053}
1054
1055CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
1056 int width )
1057{
1058 if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
1059 CountedPtr< Scantable > out = getScantable(in, false);
1060 Table& tout = out->table();
1061 out->frequencies().rescale(width, "BIN");
1062 ArrayColumn<Float> specCol(tout, "SPECTRA");
1063 ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
1064 for (uInt i=0; i < tout.nrow(); ++i ) {
1065 MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
1066 MaskedArray<Float> maout;
1067 LatticeUtilities::bin(maout, main, 0, Int(width));
1068 /// @todo implement channel based tsys binning
1069 specCol.put(i, maout.getArray());
1070 flagCol.put(i, flagsFromMA(maout));
1071 // take only the first binned spectrum's length for the deprecated
1072 // global header item nChan
1073 if (i==0) tout.rwKeywordSet().define(String("nChan"),
1074 Int(maout.getArray().nelements()));
1075 }
1076 return out;
1077}
1078
1079CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
1080 const std::string& method,
1081 float width )
1082//
1083// Should add the possibility of width being specified in km/s. This means
1084// that for each freqID (SpectralCoordinate) we will need to convert to an
1085// average channel width (say at the reference pixel). Then we would need
1086// to be careful to make sure each spectrum (of different freqID)
1087// is the same length.
1088//
1089{
1090 //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1091 Int interpMethod(stringToIMethod(method));
1092
1093 CountedPtr< Scantable > out = getScantable(in, false);
1094 Table& tout = out->table();
1095
1096// Resample SpectralCoordinates (one per freqID)
1097 out->frequencies().rescale(width, "RESAMPLE");
1098 TableIterator iter(tout, "IFNO");
1099 TableRow row(tout);
1100 while ( !iter.pastEnd() ) {
1101 Table tab = iter.table();
1102 ArrayColumn<Float> specCol(tab, "SPECTRA");
1103 //ArrayColumn<Float> tsysCol(tout, "TSYS");
1104 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1105 Vector<Float> spec;
1106 Vector<uChar> flag;
1107 specCol.get(0,spec); // the number of channels should be constant per IF
1108 uInt nChanIn = spec.nelements();
1109 Vector<Float> xIn(nChanIn); indgen(xIn);
1110 Int fac = Int(nChanIn/width);
1111 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
1112 uInt k = 0;
1113 Float x = 0.0;
1114 while (x < Float(nChanIn) ) {
1115 xOut(k) = x;
1116 k++;
1117 x += width;
1118 }
1119 uInt nChanOut = k;
1120 xOut.resize(nChanOut, True);
1121 // process all rows for this IFNO
1122 Vector<Float> specOut;
1123 Vector<Bool> maskOut;
1124 Vector<uChar> flagOut;
1125 for (uInt i=0; i < tab.nrow(); ++i) {
1126 specCol.get(i, spec);
1127 flagCol.get(i, flag);
1128 Vector<Bool> mask(flag.nelements());
1129 convertArray(mask, flag);
1130
1131 IPosition shapeIn(spec.shape());
1132 //sh.nchan = nChanOut;
1133 InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
1134 xIn, spec, mask,
1135 interpMethod, True, True);
1136 /// @todo do the same for channel based Tsys
1137 flagOut.resize(maskOut.nelements());
1138 convertArray(flagOut, maskOut);
1139 specCol.put(i, specOut);
1140 flagCol.put(i, flagOut);
1141 }
1142 ++iter;
1143 }
1144
1145 return out;
1146}
1147
1148STMath::imethod STMath::stringToIMethod(const std::string& in)
1149{
1150 static STMath::imap lookup;
1151
1152 // initialize the lookup table if necessary
1153 if ( lookup.empty() ) {
1154 lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
1155 lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
1156 lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
1157 lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
1158 }
1159
1160 STMath::imap::const_iterator iter = lookup.find(in);
1161
1162 if ( lookup.end() == iter ) {
1163 std::string message = in;
1164 message += " is not a valid interpolation mode";
1165 throw(AipsError(message));
1166 }
1167 return iter->second;
1168}
1169
1170WeightType STMath::stringToWeight(const std::string& in)
1171{
1172 static std::map<std::string, WeightType> lookup;
1173
1174 // initialize the lookup table if necessary
1175 if ( lookup.empty() ) {
1176 lookup["NONE"] = asap::NONE;
1177 lookup["TINT"] = asap::TINT;
1178 lookup["TINTSYS"] = asap::TINTSYS;
1179 lookup["TSYS"] = asap::TSYS;
1180 lookup["VAR"] = asap::VAR;
1181 }
1182
1183 std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
1184
1185 if ( lookup.end() == iter ) {
1186 std::string message = in;
1187 message += " is not a valid weighting mode";
1188 throw(AipsError(message));
1189 }
1190 return iter->second;
1191}
1192
1193CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
1194 const vector< float > & coeff,
1195 const std::string & filename,
1196 const std::string& method)
1197{
1198 // Get elevation data from Scantable and convert to degrees
1199 CountedPtr< Scantable > out = getScantable(in, false);
1200 Table& tab = out->table();
1201 ROScalarColumn<Float> elev(tab, "ELEVATION");
1202 Vector<Float> x = elev.getColumn();
1203 x *= Float(180 / C::pi); // Degrees
1204
1205 Vector<Float> coeffs(coeff);
1206 const uInt nc = coeffs.nelements();
1207 if ( filename.length() > 0 && nc > 0 ) {
1208 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
1209 }
1210
1211 // Correct
1212 if ( nc > 0 || filename.length() == 0 ) {
1213 // Find instrument
1214 Bool throwit = True;
1215 Instrument inst =
1216 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1217 throwit);
1218
1219 // Set polynomial
1220 Polynomial<Float>* ppoly = 0;
1221 Vector<Float> coeff;
1222 String msg;
1223 if ( nc > 0 ) {
1224 ppoly = new Polynomial<Float>(nc);
1225 coeff = coeffs;
1226 msg = String("user");
1227 } else {
1228 STAttr sdAttr;
1229 coeff = sdAttr.gainElevationPoly(inst);
1230 ppoly = new Polynomial<Float>(3);
1231 msg = String("built in");
1232 }
1233
1234 if ( coeff.nelements() > 0 ) {
1235 ppoly->setCoefficients(coeff);
1236 } else {
1237 delete ppoly;
1238 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
1239 }
1240 ostringstream oss;
1241 oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
1242 oss << " " << coeff;
1243 pushLog(String(oss));
1244 const uInt nrow = tab.nrow();
1245 Vector<Float> factor(nrow);
1246 for ( uInt i=0; i < nrow; ++i ) {
1247 factor[i] = 1.0 / (*ppoly)(x[i]);
1248 }
1249 delete ppoly;
1250 scaleByVector(tab, factor, true);
1251
1252 } else {
1253 // Read and correct
1254 pushLog("Making correction from ascii Table");
1255 scaleFromAsciiTable(tab, filename, method, x, true);
1256 }
1257 return out;
1258}
1259
1260void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
1261 const std::string& method,
1262 const Vector<Float>& xout, bool dotsys)
1263{
1264
1265// Read gain-elevation ascii file data into a Table.
1266
1267 String formatString;
1268 Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
1269 scaleFromTable(in, tbl, method, xout, dotsys);
1270}
1271
1272void STMath::scaleFromTable(Table& in,
1273 const Table& table,
1274 const std::string& method,
1275 const Vector<Float>& xout, bool dotsys)
1276{
1277
1278 ROScalarColumn<Float> geElCol(table, "ELEVATION");
1279 ROScalarColumn<Float> geFacCol(table, "FACTOR");
1280 Vector<Float> xin = geElCol.getColumn();
1281 Vector<Float> yin = geFacCol.getColumn();
1282 Vector<Bool> maskin(xin.nelements(),True);
1283
1284 // Interpolate (and extrapolate) with desired method
1285
1286 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1287
1288 Vector<Float> yout;
1289 Vector<Bool> maskout;
1290 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
1291 xin, yin, maskin, interp,
1292 True, True);
1293
1294 scaleByVector(in, Float(1.0)/yout, dotsys);
1295}
1296
1297void STMath::scaleByVector( Table& in,
1298 const Vector< Float >& factor,
1299 bool dotsys )
1300{
1301 uInt nrow = in.nrow();
1302 if ( factor.nelements() != nrow ) {
1303 throw(AipsError("factors.nelements() != table.nelements()"));
1304 }
1305 ArrayColumn<Float> specCol(in, "SPECTRA");
1306 ArrayColumn<uChar> flagCol(in, "FLAGTRA");
1307 ArrayColumn<Float> tsysCol(in, "TSYS");
1308 for (uInt i=0; i < nrow; ++i) {
1309 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1310 ma *= factor[i];
1311 specCol.put(i, ma.getArray());
1312 flagCol.put(i, flagsFromMA(ma));
1313 if ( dotsys ) {
1314 Vector<Float> tsys = tsysCol(i);
1315 tsys *= factor[i];
1316 tsysCol.put(i,tsys);
1317 }
1318 }
1319}
1320
1321CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
1322 float d, float etaap,
1323 float jyperk )
1324{
1325 CountedPtr< Scantable > out = getScantable(in, false);
1326 Table& tab = in->table();
1327 Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
1328 Unit K(String("K"));
1329 Unit JY(String("Jy"));
1330
1331 bool tokelvin = true;
1332 Double cfac = 1.0;
1333
1334 if ( fluxUnit == JY ) {
1335 pushLog("Converting to K");
1336 Quantum<Double> t(1.0,fluxUnit);
1337 Quantum<Double> t2 = t.get(JY);
1338 cfac = (t2 / t).getValue(); // value to Jy
1339
1340 tokelvin = true;
1341 out->setFluxUnit("K");
1342 } else if ( fluxUnit == K ) {
1343 pushLog("Converting to Jy");
1344 Quantum<Double> t(1.0,fluxUnit);
1345 Quantum<Double> t2 = t.get(K);
1346 cfac = (t2 / t).getValue(); // value to K
1347
1348 tokelvin = false;
1349 out->setFluxUnit("Jy");
1350 } else {
1351 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
1352 }
1353 // Make sure input values are converted to either Jy or K first...
1354 Float factor = cfac;
1355
1356 // Select method
1357 if (jyperk > 0.0) {
1358 factor *= jyperk;
1359 if ( tokelvin ) factor = 1.0 / jyperk;
1360 ostringstream oss;
1361 oss << "Jy/K = " << jyperk;
1362 pushLog(String(oss));
1363 Vector<Float> factors(tab.nrow(), factor);
1364 scaleByVector(tab,factors, false);
1365 } else if ( etaap > 0.0) {
1366 if (d < 0) {
1367 Instrument inst =
1368 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1369 True);
1370 STAttr sda;
1371 d = sda.diameter(inst);
1372 }
1373 jyperk = STAttr::findJyPerK(etaap, d);
1374 ostringstream oss;
1375 oss << "Jy/K = " << jyperk;
1376 pushLog(String(oss));
1377 factor *= jyperk;
1378 if ( tokelvin ) {
1379 factor = 1.0 / factor;
1380 }
1381 Vector<Float> factors(tab.nrow(), factor);
1382 scaleByVector(tab, factors, False);
1383 } else {
1384
1385 // OK now we must deal with automatic look up of values.
1386 // We must also deal with the fact that the factors need
1387 // to be computed per IF and may be different and may
1388 // change per integration.
1389
1390 pushLog("Looking up conversion factors");
1391 convertBrightnessUnits(out, tokelvin, cfac);
1392 }
1393
1394 return out;
1395}
1396
1397void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
1398 bool tokelvin, float cfac )
1399{
1400 Table& table = in->table();
1401 Instrument inst =
1402 STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
1403 TableIterator iter(table, "FREQ_ID");
1404 STFrequencies stfreqs = in->frequencies();
1405 STAttr sdAtt;
1406 while (!iter.pastEnd()) {
1407 Table tab = iter.table();
1408 ArrayColumn<Float> specCol(tab, "SPECTRA");
1409 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1410 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
1411 MEpoch::ROScalarColumn timeCol(tab, "TIME");
1412
1413 uInt freqid; freqidCol.get(0, freqid);
1414 Vector<Float> tmpspec; specCol.get(0, tmpspec);
1415 // STAttr.JyPerK has a Vector interface... change sometime.
1416 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
1417 for ( uInt i=0; i<tab.nrow(); ++i) {
1418 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
1419 Float factor = cfac * jyperk;
1420 if ( tokelvin ) factor = Float(1.0) / factor;
1421 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1422 ma *= factor;
1423 specCol.put(i, ma.getArray());
1424 flagCol.put(i, flagsFromMA(ma));
1425 }
1426 ++iter;
1427 }
1428}
1429
1430CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
1431 float tau )
1432{
1433 CountedPtr< Scantable > out = getScantable(in, false);
1434
1435 Table tab = out->table();
1436 ROScalarColumn<Float> elev(tab, "ELEVATION");
1437 ArrayColumn<Float> specCol(tab, "SPECTRA");
1438 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1439 for ( uInt i=0; i<tab.nrow(); ++i) {
1440 Float zdist = Float(C::pi_2) - elev(i);
1441 Float factor = exp(tau/cos(zdist));
1442 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1443 ma *= factor;
1444 specCol.put(i, ma.getArray());
1445 flagCol.put(i, flagsFromMA(ma));
1446 }
1447 return out;
1448}
1449
1450CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
1451 const std::string& kernel,
1452 float width )
1453{
1454 CountedPtr< Scantable > out = getScantable(in, false);
1455 Table& table = out->table();
1456 ArrayColumn<Float> specCol(table, "SPECTRA");
1457 ArrayColumn<uChar> flagCol(table, "FLAGTRA");
1458 Vector<Float> spec;
1459 Vector<uChar> flag;
1460 for ( uInt i=0; i<table.nrow(); ++i) {
1461 specCol.get(i, spec);
1462 flagCol.get(i, flag);
1463 Vector<Bool> mask(flag.nelements());
1464 convertArray(mask, flag);
1465 Vector<Float> specout;
1466 Vector<Bool> maskout;
1467 if ( kernel == "hanning" ) {
1468 mathutil::hanning(specout, maskout, spec , !mask);
1469 convertArray(flag, !maskout);
1470 } else if ( kernel == "rmedian" ) {
1471 mathutil::runningMedian(specout, maskout, spec , mask, width);
1472 convertArray(flag, maskout);
1473 }
1474 flagCol.put(i, flag);
1475 specCol.put(i, specout);
1476 }
1477 return out;
1478}
1479
1480CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
1481 const std::string& kernel, float width )
1482{
1483 if (kernel == "rmedian" || kernel == "hanning") {
1484 return smoothOther(in, kernel, width);
1485 }
1486 CountedPtr< Scantable > out = getScantable(in, false);
1487 Table& table = out->table();
1488 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
1489 // same IFNO should have same no of channels
1490 // this saves overhead
1491 TableIterator iter(table, "IFNO");
1492 while (!iter.pastEnd()) {
1493 Table tab = iter.table();
1494 ArrayColumn<Float> specCol(tab, "SPECTRA");
1495 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1496 Vector<Float> tmpspec; specCol.get(0, tmpspec);
1497 uInt nchan = tmpspec.nelements();
1498 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
1499 Convolver<Float> conv(kvec, IPosition(1,nchan));
1500 Vector<Float> spec;
1501 Vector<uChar> flag;
1502 for ( uInt i=0; i<tab.nrow(); ++i) {
1503 specCol.get(i, spec);
1504 flagCol.get(i, flag);
1505 Vector<Bool> mask(flag.nelements());
1506 convertArray(mask, flag);
1507 Vector<Float> specout;
1508 mathutil::replaceMaskByZero(specout, mask);
1509 conv.linearConv(specout, spec);
1510 specCol.put(i, specout);
1511 }
1512 ++iter;
1513 }
1514 return out;
1515}
1516
1517CountedPtr< Scantable >
1518 STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
1519{
1520 if ( in.size() < 2 ) {
1521 throw(AipsError("Need at least two scantables to perform a merge."));
1522 }
1523 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
1524 bool insitu = insitu_;
1525 setInsitu(false);
1526 CountedPtr< Scantable > out = getScantable(*it, false);
1527 setInsitu(insitu);
1528 Table& tout = out->table();
1529 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
1530 ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
1531 // Renumber SCANNO to be 0-based
1532 Vector<uInt> scannos = scannocol.getColumn();
1533 uInt offset = min(scannos);
1534 scannos -= offset;
1535 scannocol.putColumn(scannos);
1536 uInt newscanno = max(scannos)+1;
1537 ++it;
1538 while ( it != in.end() ){
1539 if ( ! (*it)->conformant(*out) ) {
1540 // non conformant.
1541 pushLog(String("Warning: Can't merge scantables as header info differs."));
1542 }
1543 out->appendToHistoryTable((*it)->history());
1544 const Table& tab = (*it)->table();
1545 TableIterator scanit(tab, "SCANNO");
1546 while (!scanit.pastEnd()) {
1547 TableIterator freqit(scanit.table(), "FREQ_ID");
1548 while ( !freqit.pastEnd() ) {
1549 Table thetab = freqit.table();
1550 uInt nrow = tout.nrow();
1551 //tout.addRow(thetab.nrow());
1552 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
1553 ROTableRow row(thetab);
1554 for ( uInt i=0; i<thetab.nrow(); ++i) {
1555 uInt k = nrow+i;
1556 scannocol.put(k, newscanno);
1557 const TableRecord& rec = row.get(i);
1558 Double rv,rp,inc;
1559 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
1560 uInt id;
1561 id = out->frequencies().addEntry(rp, rv, inc);
1562 freqidcol.put(k,id);
1563 String name,fname;Double rf;
1564 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
1565 id = out->molecules().addEntry(rf, name, fname);
1566 molidcol.put(k, id);
1567 Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
1568 (*it)->focus().getEntry(fax, ftan, frot, fhand,
1569 fmount,fuser, fxy, fxyp,
1570 rec.asuInt("FOCUS_ID"));
1571 id = out->focus().addEntry(fax, ftan, frot, fhand,
1572 fmount,fuser, fxy, fxyp);
1573 focusidcol.put(k, id);
1574 }
1575 ++freqit;
1576 }
1577 ++newscanno;
1578 ++scanit;
1579 }
1580 ++it;
1581 }
1582 return out;
1583}
1584
1585CountedPtr< Scantable >
1586 STMath::invertPhase( const CountedPtr < Scantable >& in )
1587{
1588 return applyToPol(in, &STPol::invertPhase, Float(0.0));
1589}
1590
1591CountedPtr< Scantable >
1592 STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
1593{
1594 return applyToPol(in, &STPol::rotatePhase, Float(phase));
1595}
1596
1597CountedPtr< Scantable >
1598 STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
1599{
1600 return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
1601}
1602
1603CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
1604 STPol::polOperation fptr,
1605 Float phase )
1606{
1607 CountedPtr< Scantable > out = getScantable(in, false);
1608 Table& tout = out->table();
1609 Block<String> cols(4);
1610 cols[0] = String("SCANNO");
1611 cols[1] = String("BEAMNO");
1612 cols[2] = String("IFNO");
1613 cols[3] = String("CYCLENO");
1614 TableIterator iter(tout, cols);
1615 CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
1616 out->getPolType() );
1617 while (!iter.pastEnd()) {
1618 Table t = iter.table();
1619 ArrayColumn<Float> speccol(t, "SPECTRA");
1620 ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
1621 ScalarColumn<Float> parancol(t, "PARANGLE");
1622 Matrix<Float> pols(speccol.getColumn());
1623 try {
1624 stpol->setSpectra(pols);
1625 Float fang,fhand,parang;
1626 fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
1627 fhand = in->focusTable_.getFeedHand(focidcol(0));
1628 parang = parancol(0);
1629 /// @todo re-enable this
1630 // disable total feed angle to support paralactifying Caswell style
1631 stpol->setPhaseCorrections(parang, -parang, fhand);
1632 // use a member function pointer in STPol. This only works on
1633 // the STPol pointer itself, not the Counted Pointer so
1634 // derefernce it.
1635 (&(*(stpol))->*fptr)(phase);
1636 speccol.putColumn(stpol->getSpectra());
1637 } catch (AipsError& e) {
1638 //delete stpol;stpol=0;
1639 throw(e);
1640 }
1641 ++iter;
1642 }
1643 //delete stpol;stpol=0;
1644 return out;
1645}
1646
1647CountedPtr< Scantable >
1648 STMath::swapPolarisations( const CountedPtr< Scantable > & in )
1649{
1650 CountedPtr< Scantable > out = getScantable(in, false);
1651 Table& tout = out->table();
1652 Table t0 = tout(tout.col("POLNO") == 0);
1653 Table t1 = tout(tout.col("POLNO") == 1);
1654 if ( t0.nrow() != t1.nrow() )
1655 throw(AipsError("Inconsistent number of polarisations"));
1656 ArrayColumn<Float> speccol0(t0, "SPECTRA");
1657 ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
1658 ArrayColumn<Float> speccol1(t1, "SPECTRA");
1659 ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
1660 Matrix<Float> s0 = speccol0.getColumn();
1661 Matrix<uChar> f0 = flagcol0.getColumn();
1662 speccol0.putColumn(speccol1.getColumn());
1663 flagcol0.putColumn(flagcol1.getColumn());
1664 speccol1.putColumn(s0);
1665 flagcol1.putColumn(f0);
1666 return out;
1667}
1668
1669CountedPtr< Scantable >
1670 STMath::averagePolarisations( const CountedPtr< Scantable > & in,
1671 const std::vector<bool>& mask,
1672 const std::string& weight )
1673{
1674 if (in->npol() < 2 )
1675 throw(AipsError("averagePolarisations can only be applied to two or more"
1676 "polarisations"));
1677 bool insitu = insitu_;
1678 setInsitu(false);
1679 CountedPtr< Scantable > pols = getScantable(in, true);
1680 setInsitu(insitu);
1681 Table& tout = pols->table();
1682 std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
1683 Table tab = tableCommand(taql, in->table());
1684 if (tab.nrow() == 0 )
1685 throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
1686 TableCopy::copyRows(tout, tab);
1687 TableVector<uInt> vec(tout, "POLNO");
1688 vec = 0;
1689 pols->table_.rwKeywordSet().define("nPol", Int(1));
1690 //pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
1691 pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
1692 std::vector<CountedPtr<Scantable> > vpols;
1693 vpols.push_back(pols);
1694 CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
1695 return out;
1696}
1697
1698CountedPtr< Scantable >
1699 STMath::averageBeams( const CountedPtr< Scantable > & in,
1700 const std::vector<bool>& mask,
1701 const std::string& weight )
1702{
1703 bool insitu = insitu_;
1704 setInsitu(false);
1705 CountedPtr< Scantable > beams = getScantable(in, false);
1706 setInsitu(insitu);
1707 Table& tout = beams->table();
1708 // give all rows the same BEAMNO
1709 TableVector<uInt> vec(tout, "BEAMNO");
1710 vec = 0;
1711 beams->table_.rwKeywordSet().define("nBeam", Int(1));
1712 std::vector<CountedPtr<Scantable> > vbeams;
1713 vbeams.push_back(beams);
1714 CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
1715 return out;
1716}
1717
1718
1719CountedPtr< Scantable >
1720 asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
1721 const std::string & refTime,
1722 const std::string & method)
1723{
1724 // clone as this is not working insitu
1725 bool insitu = insitu_;
1726 setInsitu(false);
1727 CountedPtr< Scantable > out = getScantable(in, false);
1728 setInsitu(insitu);
1729 Table& tout = out->table();
1730 // Get reference Epoch to time of first row or given String
1731 Unit DAY(String("d"));
1732 MEpoch::Ref epochRef(in->getTimeReference());
1733 MEpoch refEpoch;
1734 if (refTime.length()>0) {
1735 Quantum<Double> qt;
1736 if (MVTime::read(qt,refTime)) {
1737 MVEpoch mv(qt);
1738 refEpoch = MEpoch(mv, epochRef);
1739 } else {
1740 throw(AipsError("Invalid format for Epoch string"));
1741 }
1742 } else {
1743 refEpoch = in->timeCol_(0);
1744 }
1745 MPosition refPos = in->getAntennaPosition();
1746
1747 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1748 // test if user frame is different to base frame
1749 if ( in->frequencies().getFrameString(true)
1750 == in->frequencies().getFrameString(false) ) {
1751 throw(AipsError("Can't convert as no output frame has been set"
1752 " (use set_freqframe) or it is aligned already."));
1753 }
1754 MFrequency::Types system = in->frequencies().getFrame();
1755 MVTime mvt(refEpoch.getValue());
1756 String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
1757 ostringstream oss;
1758 oss << "Aligned at reference Epoch " << epochout
1759 << " in frame " << MFrequency::showType(system);
1760 pushLog(String(oss));
1761 // set up the iterator
1762 Block<String> cols(4);
1763 // select by constant direction
1764 cols[0] = String("SRCNAME");
1765 cols[1] = String("BEAMNO");
1766 // select by IF ( no of channels varies over this )
1767 cols[2] = String("IFNO");
1768 // select by restfrequency
1769 cols[3] = String("MOLECULE_ID");
1770 TableIterator iter(tout, cols);
1771 while ( !iter.pastEnd() ) {
1772 Table t = iter.table();
1773 MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1774 TableIterator fiter(t, "FREQ_ID");
1775 // determine nchan from the first row. This should work as
1776 // we are iterating over BEAMNO and IFNO // we should have constant direction
1777
1778 ROArrayColumn<Float> sCol(t, "SPECTRA");
1779 MDirection direction = dirCol(0);
1780 uInt nchan = sCol(0).nelements();
1781 while ( !fiter.pastEnd() ) {
1782 Table ftab = fiter.table();
1783 ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
1784 // get the SpectralCoordinate for the freqid, which we are iterating over
1785 SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1786 FrequencyAligner<Float> fa( sC, nchan, refEpoch,
1787 direction, refPos, system );
1788 // realign the SpectralCoordinate and put into the output Scantable
1789 Vector<String> units(1);
1790 units = String("Hz");
1791 Bool linear=True;
1792 SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1793 sc2.setWorldAxisUnits(units);
1794 uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
1795 sc2.referenceValue()[0],
1796 sc2.increment()[0]);
1797 TableVector<uInt> tvec(ftab, "FREQ_ID");
1798 tvec = id;
1799 // create the "global" abcissa for alignment with same FREQ_ID
1800 Vector<Double> abc(nchan);
1801 Double w;
1802 for (uInt i=0; i<nchan; i++) {
1803 sC.toWorld(w,Double(i));
1804 abc[i] = w;
1805 }
1806 // cache abcissa for same time stamps, so iterate over those
1807 TableIterator timeiter(ftab, "TIME");
1808 while ( !timeiter.pastEnd() ) {
1809 Table tab = timeiter.table();
1810 ArrayColumn<Float> specCol(tab, "SPECTRA");
1811 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1812 MEpoch::ROScalarColumn timeCol(tab, "TIME");
1813 // use align abcissa cache after the first row
1814 bool first = true;
1815 // these rows should be just be POLNO
1816 for (int i=0; i<int(tab.nrow()); ++i) {
1817 // input values
1818 Vector<uChar> flag = flagCol(i);
1819 Vector<Bool> mask(flag.shape());
1820 Vector<Float> specOut, spec;
1821 spec = specCol(i);
1822 Vector<Bool> maskOut;Vector<uChar> flagOut;
1823 convertArray(mask, flag);
1824 // alignment
1825 Bool ok = fa.align(specOut, maskOut, abc, spec,
1826 mask, timeCol(i), !first,
1827 interp, False);
1828 // back into scantable
1829 flagOut.resize(maskOut.nelements());
1830 convertArray(flagOut, maskOut);
1831 flagCol.put(i, flagOut);
1832 specCol.put(i, specOut);
1833 // start abcissa caching
1834 first = false;
1835 }
1836 // next timestamp
1837 ++timeiter;
1838 }
1839 // next FREQ_ID
1840 ++fiter;
1841 }
1842 // next aligner
1843 ++iter;
1844 }
1845 // set this afterwards to ensure we are doing insitu correctly.
1846 out->frequencies().setFrame(system, true);
1847 return out;
1848}
1849
1850CountedPtr<Scantable>
1851 asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
1852 const std::string & newtype )
1853{
1854 if (in->npol() != 2 && in->npol() != 4)
1855 throw(AipsError("Can only convert two or four polarisations."));
1856 if ( in->getPolType() == newtype )
1857 throw(AipsError("No need to convert."));
1858 if ( ! in->selector_.empty() )
1859 throw(AipsError("Can only convert whole scantable. Unset the selection."));
1860 bool insitu = insitu_;
1861 setInsitu(false);
1862 CountedPtr< Scantable > out = getScantable(in, true);
1863 setInsitu(insitu);
1864 Table& tout = out->table();
1865 tout.rwKeywordSet().define("POLTYPE", String(newtype));
1866
1867 Block<String> cols(4);
1868 cols[0] = "SCANNO";
1869 cols[1] = "CYCLENO";
1870 cols[2] = "BEAMNO";
1871 cols[3] = "IFNO";
1872 TableIterator it(in->originalTable_, cols);
1873 String basetype = in->getPolType();
1874 STPol* stpol = STPol::getPolClass(in->factories_, basetype);
1875 try {
1876 while ( !it.pastEnd() ) {
1877 Table tab = it.table();
1878 uInt row = tab.rowNumbers()[0];
1879 stpol->setSpectra(in->getPolMatrix(row));
1880 Float fang,fhand,parang;
1881 fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
1882 fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
1883 parang = in->paraCol_(row);
1884 /// @todo re-enable this
1885 // disable total feed angle to support paralactifying Caswell style
1886 stpol->setPhaseCorrections(parang, -parang, fhand);
1887 Int npolout = 0;
1888 for (uInt i=0; i<tab.nrow(); ++i) {
1889 Vector<Float> outvec = stpol->getSpectrum(i, newtype);
1890 if ( outvec.nelements() > 0 ) {
1891 tout.addRow();
1892 TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
1893 ArrayColumn<Float> sCol(tout,"SPECTRA");
1894 ScalarColumn<uInt> pCol(tout,"POLNO");
1895 sCol.put(tout.nrow()-1 ,outvec);
1896 pCol.put(tout.nrow()-1 ,uInt(npolout));
1897 npolout++;
1898 }
1899 }
1900 tout.rwKeywordSet().define("nPol", npolout);
1901 ++it;
1902 }
1903 } catch (AipsError& e) {
1904 delete stpol;
1905 throw(e);
1906 }
1907 delete stpol;
1908 return out;
1909}
1910
1911CountedPtr< Scantable >
1912 asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
1913 const std::string & scantype )
1914{
1915 bool insitu = insitu_;
1916 setInsitu(false);
1917 CountedPtr< Scantable > out = getScantable(in, true);
1918 setInsitu(insitu);
1919 Table& tout = out->table();
1920 std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
1921 if (scantype == "on") {
1922 taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
1923 }
1924 Table tab = tableCommand(taql, in->table());
1925 TableCopy::copyRows(tout, tab);
1926 if (scantype == "on") {
1927 // re-index SCANNO to 0
1928 TableVector<uInt> vec(tout, "SCANNO");
1929 vec = 0;
1930 }
1931 return out;
1932}
1933
1934CountedPtr< Scantable >
1935 asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
1936 double frequency, double width )
1937{
1938 CountedPtr< Scantable > out = getScantable(in, false);
1939 Table& tout = out->table();
1940 TableIterator iter(tout, "FREQ_ID");
1941 FFTServer<Float,Complex> ffts;
1942 while ( !iter.pastEnd() ) {
1943 Table tab = iter.table();
1944 Double rp,rv,inc;
1945 ROTableRow row(tab);
1946 const TableRecord& rec = row.get(0);
1947 uInt freqid = rec.asuInt("FREQ_ID");
1948 out->frequencies().getEntry(rp, rv, inc, freqid);
1949 ArrayColumn<Float> specCol(tab, "SPECTRA");
1950 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1951 for (int i=0; i<int(tab.nrow()); ++i) {
1952 Vector<Float> spec = specCol(i);
1953 Vector<uChar> flag = flagCol(i);
1954 Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
1955 Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
1956 for (int k=0; k < flag.nelements(); ++k ) {
1957 if (flag[k] > 0) {
1958 spec[k] = 0.0;
1959 }
1960 }
1961 Vector<Complex> lags;
1962 ffts.fft0(lags, spec);
1963 Int start = max(0, lag0);
1964 Int end = min(Int(lags.nelements()-1), lag1);
1965 if (start == end) {
1966 lags[start] = Complex(0.0);
1967 } else {
1968 for (int j=start; j <=end ;++j) {
1969 lags[j] = Complex(0.0);
1970 }
1971 }
1972 ffts.fft0(spec, lags);
1973 specCol.put(i, spec);
1974 }
1975 ++iter;
1976 }
1977 return out;
1978}
Note: See TracBrowser for help on using the repository browser.