source: trunk/src/STMath.cpp@ 1412

Last change on this file since 1412 was 1391, checked in by Malte Marquarding, 17 years ago

merge from alma branch to get alma/GBT support. Commented out fluxUnit changes as they are using a chnaged interface to PKSreader/writer. Also commented out progress meter related code.

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