source: trunk/src/STMath.cpp@ 1817

Last change on this file since 1817 was 1689, checked in by Malte Marquarding, 15 years ago

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