source: trunk/src/SDMath.cc@ 528

Last change on this file since 528 was 519, checked in by kil064, 20 years ago

fix typo
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 58.9 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDMath.cc: A collection of single dish mathematical operations
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id:
30//#---------------------------------------------------------------------------
31#include <vector>
32
33#include <casa/aips.h>
34#include <casa/iostream.h>
35#include <casa/iomanip.h>
36#include <casa/BasicSL/String.h>
37#include <casa/Arrays/IPosition.h>
38#include <casa/Arrays/Array.h>
39#include <casa/Arrays/ArrayIter.h>
40#include <casa/Arrays/VectorIter.h>
41#include <casa/Arrays/ArrayMath.h>
42#include <casa/Arrays/ArrayLogical.h>
43#include <casa/Arrays/MaskedArray.h>
44#include <casa/Arrays/MaskArrMath.h>
45#include <casa/Arrays/MaskArrLogi.h>
46#include <casa/Arrays/Matrix.h>
47#include <casa/BasicMath/Math.h>
48#include <casa/Containers/Block.h>
49#include <casa/Exceptions.h>
50#include <casa/Quanta/Quantum.h>
51#include <casa/Quanta/Unit.h>
52#include <casa/Quanta/MVEpoch.h>
53#include <casa/Quanta/MVTime.h>
54#include <casa/Utilities/Assert.h>
55
56#include <coordinates/Coordinates/SpectralCoordinate.h>
57#include <coordinates/Coordinates/CoordinateSystem.h>
58#include <coordinates/Coordinates/CoordinateUtil.h>
59#include <coordinates/Coordinates/FrequencyAligner.h>
60
61#include <lattices/Lattices/LatticeUtilities.h>
62#include <lattices/Lattices/RebinLattice.h>
63
64#include <measures/Measures/MEpoch.h>
65#include <measures/Measures/MDirection.h>
66#include <measures/Measures/MPosition.h>
67
68#include <scimath/Mathematics/VectorKernel.h>
69#include <scimath/Mathematics/Convolver.h>
70#include <scimath/Mathematics/InterpolateArray1D.h>
71#include <scimath/Functionals/Polynomial.h>
72
73#include <tables/Tables/Table.h>
74#include <tables/Tables/ScalarColumn.h>
75#include <tables/Tables/ArrayColumn.h>
76#include <tables/Tables/ReadAsciiTable.h>
77
78#include "MathUtils.h"
79#include "SDDefs.h"
80#include "SDAttr.h"
81#include "SDContainer.h"
82#include "SDMemTable.h"
83
84#include "SDMath.h"
85#include "SDPol.h"
86
87using namespace casa;
88using namespace asap;
89
90
91SDMath::SDMath()
92{;}
93
94SDMath::SDMath(const SDMath& other)
95{
96
97// No state
98
99}
100
101SDMath& SDMath::operator=(const SDMath& other)
102{
103 if (this != &other) {
104// No state
105 }
106 return *this;
107}
108
109SDMath::~SDMath()
110{;}
111
112
113
114SDMemTable* SDMath::frequencyAlignment(const SDMemTable& in,
115 const String& refTime,
116 const String& method,
117 Bool perFreqID) const
118{
119// Get frame info from Table
120
121 std::vector<std::string> info = in.getCoordInfo();
122
123// Parse frequency system
124
125 String systemStr(info[1]);
126 String baseSystemStr(info[3]);
127 if (baseSystemStr==systemStr) {
128 throw(AipsError("You have not set a frequency frame different from the initial - use function set_freqframe"));
129 }
130//
131 MFrequency::Types freqSystem;
132 MFrequency::getType(freqSystem, systemStr);
133
134// Do it
135
136 return frequencyAlign(in, freqSystem, refTime, method, perFreqID);
137}
138
139
140
141CountedPtr<SDMemTable> SDMath::average(const Block<CountedPtr<SDMemTable> >& in,
142 const Vector<Bool>& mask, Bool scanAv,
143 const String& weightStr, Bool alignFreq) const
144//
145// Weighted averaging of spectra from one or more Tables.
146//
147{
148
149// Convert weight type
150
151 WeightType wtType = NONE;
152 convertWeightString(wtType, weightStr, True);
153
154// Create output Table by cloning from the first table
155
156 SDMemTable* pTabOut = new SDMemTable(*in[0],True);
157 if (in.nelements() > 1) {
158 for (uInt i=1; i < in.nelements(); ++i) {
159 pTabOut->appendToHistoryTable(in[i]->getHistoryTable());
160 }
161 }
162// Setup
163
164 IPosition shp = in[0]->rowAsMaskedArray(0).shape(); // Must not change
165 Array<Float> arr(shp);
166 Array<Bool> barr(shp);
167 const Bool useMask = (mask.nelements() == shp(asap::ChanAxis));
168
169// Columns from Tables
170
171 ROArrayColumn<Float> tSysCol;
172 ROScalarColumn<Double> mjdCol;
173 ROScalarColumn<String> srcNameCol;
174 ROScalarColumn<Double> intCol;
175 ROArrayColumn<uInt> fqIDCol;
176 ROScalarColumn<Int> scanIDCol;
177
178// Create accumulation MaskedArray. We accumulate for each channel,if,pol,beam
179// Note that the mask of the accumulation array will ALWAYS remain ALL True.
180// The MA is only used so that when data which is masked Bad is added to it,
181// that data does not contribute.
182
183 Array<Float> zero(shp);
184 zero=0.0;
185 Array<Bool> good(shp);
186 good = True;
187 MaskedArray<Float> sum(zero,good);
188
189// Counter arrays
190
191 Array<Float> nPts(shp); // Number of points
192 nPts = 0.0;
193 Array<Float> nInc(shp); // Increment
194 nInc = 1.0;
195
196// Create accumulation Array for variance. We accumulate for
197// each if,pol,beam, but average over channel. So we need
198// a shape with one less axis dropping channels.
199
200 const uInt nAxesSub = shp.nelements() - 1;
201 IPosition shp2(nAxesSub);
202 for (uInt i=0,j=0; i<(nAxesSub+1); i++) {
203 if (i!=asap::ChanAxis) {
204 shp2(j) = shp(i);
205 j++;
206 }
207 }
208 Array<Float> sumSq(shp2);
209 sumSq = 0.0;
210 IPosition pos2(nAxesSub,0); // For indexing
211
212// Time-related accumulators
213
214 Double time;
215 Double timeSum = 0.0;
216 Double intSum = 0.0;
217 Double interval = 0.0;
218
219// To get the right shape for the Tsys accumulator we need to
220// access a column from the first table. The shape of this
221// array must not change. Note however that since the TSysSqSum
222// array is used in a normalization process, and that I ignore the
223// channel axis replication of values for now, it loses a dimension
224
225 Array<Float> tSysSum, tSysSqSum;
226 {
227 const Table& tabIn = in[0]->table();
228 tSysCol.attach(tabIn,"TSYS");
229 tSysSum.resize(tSysCol.shape(0));
230//
231 tSysSqSum.resize(shp2);
232 }
233 tSysSum =0.0;
234 tSysSqSum = 0.0;
235 Array<Float> tSys;
236
237// Scan and row tracking
238
239 Int oldScanID = 0;
240 Int outScanID = 0;
241 Int scanID = 0;
242 Int rowStart = 0;
243 Int nAccum = 0;
244 Int tableStart = 0;
245
246// Source and FreqID
247
248 String sourceName, oldSourceName, sourceNameStart;
249 Vector<uInt> freqID, freqIDStart, oldFreqID;
250
251// Loop over tables
252
253 Float fac = 1.0;
254 const uInt nTables = in.nelements();
255 for (uInt iTab=0; iTab<nTables; iTab++) {
256
257// Should check that the frequency tables don't change if doing FreqAlignment
258
259// Attach columns to Table
260
261 const Table& tabIn = in[iTab]->table();
262 tSysCol.attach(tabIn, "TSYS");
263 mjdCol.attach(tabIn, "TIME");
264 srcNameCol.attach(tabIn, "SRCNAME");
265 intCol.attach(tabIn, "INTERVAL");
266 fqIDCol.attach(tabIn, "FREQID");
267 scanIDCol.attach(tabIn, "SCANID");
268
269// Find list of start/end rows for each scan
270
271// Loop over rows in Table
272
273 const uInt nRows = in[iTab]->nRow();
274 for (uInt iRow=0; iRow<nRows; iRow++) {
275
276// Check conformance
277
278 IPosition shp2 = in[iTab]->rowAsMaskedArray(iRow).shape();
279 if (!shp.isEqual(shp2)) {
280 throw (AipsError("Shapes for all rows must be the same"));
281 }
282
283// If we are not doing scan averages, make checks for source and
284// frequency setup and warn if averaging across them
285
286 scanIDCol.getScalar(iRow, scanID);
287
288// Get quantities from columns
289
290 srcNameCol.getScalar(iRow, sourceName);
291 mjdCol.get(iRow, time);
292 tSysCol.get(iRow, tSys);
293 intCol.get(iRow, interval);
294 fqIDCol.get(iRow, freqID);
295
296// Initialize first source and freqID
297
298 if (iRow==0 && iTab==0) {
299 sourceNameStart = sourceName;
300 freqIDStart = freqID;
301 }
302
303// If we are doing scan averages, see if we are at the end of an
304// accumulation period (scan). We must check soutce names too,
305// since we might have two tables with one scan each but different
306// source names; we shouldn't average different sources together
307
308 if (scanAv && ( (scanID != oldScanID) ||
309 (iRow==0 && iTab>0 && sourceName!=oldSourceName))) {
310
311// Normalize data in 'sum' accumulation array according to weighting scheme
312
313 normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, asap::ChanAxis, nAxesSub);
314
315// Get ScanContainer for the first row of this averaged Scan
316
317 SDContainer scOut = in[iTab]->getSDContainer(rowStart);
318
319// Fill scan container. The source and freqID come from the
320// first row of the first table that went into this average (
321// should be the same for all rows in the scan average)
322
323 Float nR(nAccum);
324 fillSDC(scOut, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
325 timeSum/nR, intSum, sourceNameStart, freqIDStart);
326
327// Write container out to Table
328
329 pTabOut->putSDContainer(scOut);
330
331// Reset accumulators
332
333 sum = 0.0;
334 sumSq = 0.0;
335 nAccum = 0;
336//
337 tSysSum =0.0;
338 tSysSqSum =0.0;
339 timeSum = 0.0;
340 intSum = 0.0;
341 nPts = 0.0;
342
343// Increment
344
345 rowStart = iRow; // First row for next accumulation
346 tableStart = iTab; // First table for next accumulation
347 sourceNameStart = sourceName; // First source name for next accumulation
348 freqIDStart = freqID; // First FreqID for next accumulation
349//
350 oldScanID = scanID;
351 outScanID += 1; // Scan ID for next accumulation period
352 }
353
354// Accumulate
355
356 accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts, tSysSum, tSysSqSum,
357 tSys, nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis,
358 nAxesSub, useMask, wtType);
359//
360 oldSourceName = sourceName;
361 oldFreqID = freqID;
362 }
363 }
364
365// OK at this point we have accumulation data which is either
366// - accumulated from all tables into one row
367// or
368// - accumulated from the last scan average
369//
370// Normalize data in 'sum' accumulation array according to weighting scheme
371
372 normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, asap::ChanAxis, nAxesSub);
373
374// Create and fill container. The container we clone will be from
375// the last Table and the first row that went into the current
376// accumulation. It probably doesn't matter that much really...
377
378 Float nR(nAccum);
379 SDContainer scOut = in[tableStart]->getSDContainer(rowStart);
380 fillSDC(scOut, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
381 timeSum/nR, intSum, sourceNameStart, freqIDStart);
382 pTabOut->putSDContainer(scOut);
383 pTabOut->resetCursor();
384//
385 return CountedPtr<SDMemTable>(pTabOut);
386}
387
388
389
390CountedPtr<SDMemTable> SDMath::binaryOperate(const CountedPtr<SDMemTable>&
391 left,
392 const CountedPtr<SDMemTable>&
393 right,
394 const String& op, Bool preserve,
395 Bool doTSys) const
396{
397
398// Check operator
399
400 String op2(op);
401 op2.upcase();
402 uInt what = 0;
403 if (op2=="ADD") {
404 what = 0;
405 } else if (op2=="SUB") {
406 what = 1;
407 } else if (op2=="MUL") {
408 what = 2;
409 } else if (op2=="DIV") {
410 what = 3;
411 } else if (op2=="QUOTIENT") {
412 what = 4;
413 doTSys = True;
414 } else {
415 throw( AipsError("Unrecognized operation"));
416 }
417
418// Check rows
419
420 const uInt nRowLeft = left->nRow();
421 const uInt nRowRight = right->nRow();
422 Bool ok = (nRowRight==1&&nRowLeft>0) ||
423 (nRowRight>=1&&nRowLeft==nRowRight);
424 if (!ok) {
425 throw (AipsError("The right Scan Table can have one row or the same number of rows as the left Scan Table"));
426 }
427
428// Input Tables
429
430 const Table& tLeft = left->table();
431 const Table& tRight = right->table();
432
433// TSys columns
434
435 ROArrayColumn<Float> tSysLeftCol, tSysRightCol;
436 if (doTSys) {
437 tSysLeftCol.attach(tLeft, "TSYS");
438 tSysRightCol.attach(tRight, "TSYS");
439 }
440
441// First row for right
442
443 Array<Float> tSysLeftArr, tSysRightArr;
444 if (doTSys) tSysRightCol.get(0, tSysRightArr);
445 MaskedArray<Float>* pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(0));
446 IPosition shpRight = pMRight->shape();
447
448// Output Table cloned from left
449
450 SDMemTable* pTabOut = new SDMemTable(*left, True);
451 pTabOut->appendToHistoryTable(right->getHistoryTable());
452// Loop over rows
453
454 for (uInt i=0; i<nRowLeft; i++) {
455
456// Get data
457
458 MaskedArray<Float> mLeft(left->rowAsMaskedArray(i));
459 IPosition shpLeft = mLeft.shape();
460 if (doTSys) tSysLeftCol.get(i, tSysLeftArr);
461//
462 if (nRowRight>1) {
463 delete pMRight;
464 pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(i));
465 shpRight = pMRight->shape();
466 if (doTSys) tSysRightCol.get(i, tSysRightArr);
467 }
468//
469 if (!shpRight.isEqual(shpLeft)) {
470 throw(AipsError("left and right scan tables are not conformant"));
471 }
472 if (doTSys) {
473 if (!tSysRightArr.shape().isEqual(tSysRightArr.shape())) {
474 throw(AipsError("left and right Tsys data are not conformant"));
475 }
476 if (!shpRight.isEqual(tSysRightArr.shape())) {
477 throw(AipsError("left and right scan tables are not conformant"));
478 }
479 }
480
481// Make container
482
483 SDContainer sc = left->getSDContainer(i);
484
485// Operate on data and TSys
486
487 if (what==0) {
488 MaskedArray<Float> tmp = mLeft + *pMRight;
489 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
490 if (doTSys) sc.putTsys(tSysLeftArr+tSysRightArr);
491 } else if (what==1) {
492 MaskedArray<Float> tmp = mLeft - *pMRight;
493 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
494 if (doTSys) sc.putTsys(tSysLeftArr-tSysRightArr);
495 } else if (what==2) {
496 MaskedArray<Float> tmp = mLeft * *pMRight;
497 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
498 if (doTSys) sc.putTsys(tSysLeftArr*tSysRightArr);
499 } else if (what==3) {
500 MaskedArray<Float> tmp = mLeft / *pMRight;
501 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
502 if (doTSys) sc.putTsys(tSysLeftArr/tSysRightArr);
503 } else if (what==4) {
504 if (preserve) {
505 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) -
506 tSysRightArr;
507 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
508 } else {
509 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) -
510 tSysLeftArr;
511 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
512 }
513 sc.putTsys(tSysRightArr);
514 }
515
516// Put new row in output Table
517
518 pTabOut->putSDContainer(sc);
519 }
520 if (pMRight) delete pMRight;
521 pTabOut->resetCursor();
522
523 return CountedPtr<SDMemTable>(pTabOut);
524}
525
526
527
528std::vector<float> SDMath::statistic(const CountedPtr<SDMemTable>& in,
529 const Vector<Bool>& mask,
530 const String& which, Int row) const
531//
532// Perhaps iteration over pol/beam/if should be in here
533// and inside the nrow iteration ?
534//
535{
536 const uInt nRow = in->nRow();
537
538// Specify cursor location
539
540 IPosition start, end;
541 Bool doAll = False;
542 setCursorSlice (start, end, doAll, *in);
543
544// Loop over rows
545
546 const uInt nEl = mask.nelements();
547 uInt iStart = 0;
548 uInt iEnd = in->nRow()-1;
549//
550 if (row>=0) {
551 iStart = row;
552 iEnd = row;
553 }
554//
555 std::vector<float> result(iEnd-iStart+1);
556 for (uInt ii=iStart; ii <= iEnd; ++ii) {
557
558// Get row and deconstruct
559
560 MaskedArray<Float> dataIn = (in->rowAsMaskedArray(ii))(start,end);
561 Array<Float> v = dataIn.getArray().nonDegenerate();
562 Array<Bool> m = dataIn.getMask().nonDegenerate();
563
564// Access desired piece of data
565
566// Array<Float> v((arr(start,end)).nonDegenerate());
567// Array<Bool> m((barr(start,end)).nonDegenerate());
568
569// Apply OTF mask
570
571 MaskedArray<Float> tmp;
572 if (m.nelements()==nEl) {
573 tmp.setData(v,m&&mask);
574 } else {
575 tmp.setData(v,m);
576 }
577
578// Get statistic
579
580 result[ii-iStart] = mathutil::statistics(which, tmp);
581 }
582//
583 return result;
584}
585
586
587SDMemTable* SDMath::bin(const SDMemTable& in, Int width) const
588{
589 SDHeader sh = in.getSDHeader();
590 SDMemTable* pTabOut = new SDMemTable(in, True);
591
592// Bin up SpectralCoordinates
593
594 IPosition factors(1);
595 factors(0) = width;
596 for (uInt j=0; j<in.nCoordinates(); ++j) {
597 CoordinateSystem cSys;
598 cSys.addCoordinate(in.getSpectralCoordinate(j));
599 CoordinateSystem cSysBin =
600 CoordinateUtil::makeBinnedCoordinateSystem(factors, cSys, False);
601//
602 SpectralCoordinate sCBin = cSysBin.spectralCoordinate(0);
603 pTabOut->setCoordinate(sCBin, j);
604 }
605
606// Use RebinLattice to find shape
607
608 IPosition shapeIn(1,sh.nchan);
609 IPosition shapeOut = RebinLattice<Float>::rebinShape(shapeIn, factors);
610 sh.nchan = shapeOut(0);
611 pTabOut->putSDHeader(sh);
612
613// Loop over rows and bin along channel axis
614
615 for (uInt i=0; i < in.nRow(); ++i) {
616 SDContainer sc = in.getSDContainer(i);
617//
618 Array<Float> tSys(sc.getTsys()); // Get it out before sc changes shape
619
620// Bin up spectrum
621
622 MaskedArray<Float> marr(in.rowAsMaskedArray(i));
623 MaskedArray<Float> marrout;
624 LatticeUtilities::bin(marrout, marr, asap::ChanAxis, width);
625
626// Put back the binned data and flags
627
628 IPosition ip2 = marrout.shape();
629 sc.resize(ip2);
630//
631 putDataInSDC(sc, marrout.getArray(), marrout.getMask());
632
633// Bin up Tsys.
634
635 Array<Bool> allGood(tSys.shape(),True);
636 MaskedArray<Float> tSysIn(tSys, allGood, True);
637//
638 MaskedArray<Float> tSysOut;
639 LatticeUtilities::bin(tSysOut, tSysIn, asap::ChanAxis, width);
640 sc.putTsys(tSysOut.getArray());
641//
642 pTabOut->putSDContainer(sc);
643 }
644 return pTabOut;
645}
646
647SDMemTable* SDMath::resample(const SDMemTable& in, const String& methodStr,
648 Float width) const
649//
650// Should add the possibility of width being specified in km/s. This means
651// that for each freqID (SpectralCoordinate) we will need to convert to an
652// average channel width (say at the reference pixel). Then we would need
653// to be careful to make sure each spectrum (of different freqID)
654// is the same length.
655//
656{
657 Bool doVel = False;
658 if (doVel) {
659 for (uInt j=0; j<in.nCoordinates(); ++j) {
660 SpectralCoordinate sC = in.getSpectralCoordinate(j);
661 }
662 }
663
664// Interpolation method
665
666 InterpolateArray1D<Double,Float>::InterpolationMethod interp;
667 convertInterpString(interp, methodStr);
668 Int interpMethod(interp);
669
670// Make output table
671
672 SDMemTable* pTabOut = new SDMemTable(in, True);
673
674// Resample SpectralCoordinates (one per freqID)
675
676 const uInt nCoord = in.nCoordinates();
677 Vector<Float> offset(1,0.0);
678 Vector<Float> factors(1,1.0/width);
679 Vector<Int> newShape;
680 for (uInt j=0; j<in.nCoordinates(); ++j) {
681 CoordinateSystem cSys;
682 cSys.addCoordinate(in.getSpectralCoordinate(j));
683 CoordinateSystem cSys2 = cSys.subImage(offset, factors, newShape);
684 SpectralCoordinate sC = cSys2.spectralCoordinate(0);
685//
686 pTabOut->setCoordinate(sC, j);
687 }
688
689// Get header
690
691 SDHeader sh = in.getSDHeader();
692
693// Generate resampling vectors
694
695 const uInt nChanIn = sh.nchan;
696 Vector<Float> xIn(nChanIn);
697 indgen(xIn);
698//
699 Int fac = Int(nChanIn/width);
700 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
701 uInt i = 0;
702 Float x = 0.0;
703 Bool more = True;
704 while (more) {
705 xOut(i) = x;
706//
707 i++;
708 x += width;
709 if (x>nChanIn-1) more = False;
710 }
711 const uInt nChanOut = i;
712 xOut.resize(nChanOut,True);
713//
714 IPosition shapeIn(in.rowAsMaskedArray(0).shape());
715 sh.nchan = nChanOut;
716 pTabOut->putSDHeader(sh);
717
718// Loop over rows and resample along channel axis
719
720 Array<Float> valuesOut;
721 Array<Bool> maskOut;
722 Array<Float> tSysOut;
723 Array<Bool> tSysMaskIn(shapeIn,True);
724 Array<Bool> tSysMaskOut;
725 for (uInt i=0; i < in.nRow(); ++i) {
726
727// Get container
728
729 SDContainer sc = in.getSDContainer(i);
730
731// Get data and Tsys
732
733 const Array<Float>& tSysIn = sc.getTsys();
734 const MaskedArray<Float>& dataIn(in.rowAsMaskedArray(i));
735 Array<Float> valuesIn = dataIn.getArray();
736 Array<Bool> maskIn = dataIn.getMask();
737
738// Interpolate data
739
740 InterpolateArray1D<Float,Float>::interpolate(valuesOut, maskOut, xOut,
741 xIn, valuesIn, maskIn,
742 interpMethod, True, True);
743 sc.resize(valuesOut.shape());
744 putDataInSDC(sc, valuesOut, maskOut);
745
746// Interpolate TSys
747
748 InterpolateArray1D<Float,Float>::interpolate(tSysOut, tSysMaskOut, xOut,
749 xIn, tSysIn, tSysMaskIn,
750 interpMethod, True, True);
751 sc.putTsys(tSysOut);
752
753// Put container in output
754
755 pTabOut->putSDContainer(sc);
756 }
757//
758 return pTabOut;
759}
760
761SDMemTable* SDMath::unaryOperate(const SDMemTable& in, Float val, Bool doAll,
762 uInt what, Bool doTSys) const
763//
764// what = 0 Multiply
765// 1 Add
766{
767 SDMemTable* pOut = new SDMemTable(in,False);
768 const Table& tOut = pOut->table();
769 ArrayColumn<Float> specCol(tOut,"SPECTRA");
770 ArrayColumn<Float> tSysCol(tOut,"TSYS");
771 Array<Float> tSysArr;
772
773// Get data slice bounds
774
775 IPosition start, end;
776 setCursorSlice (start, end, doAll, in);
777//
778 for (uInt i=0; i<tOut.nrow(); i++) {
779
780// Modify data
781
782 MaskedArray<Float> dataIn(pOut->rowAsMaskedArray(i));
783 MaskedArray<Float> dataIn2 = dataIn(start,end); // Reference
784 if (what==0) {
785 dataIn2 *= val;
786 } else if (what==1) {
787 dataIn2 += val;
788 }
789 specCol.put(i, dataIn.getArray());
790
791// Modify Tsys
792
793 if (doTSys) {
794 tSysCol.get(i, tSysArr);
795 Array<Float> tSysArr2 = tSysArr(start,end); // Reference
796 if (what==0) {
797 tSysArr2 *= val;
798 } else if (what==1) {
799 tSysArr2 += val;
800 }
801 tSysCol.put(i, tSysArr);
802 }
803 }
804//
805 return pOut;
806}
807
808SDMemTable* SDMath::averagePol(const SDMemTable& in, const Vector<Bool>& mask,
809 const String& weightStr) const
810//
811// Average all polarizations together, weighted by variance
812//
813{
814 WeightType wtType = NONE;
815 convertWeightString(wtType, weightStr, False);
816
817 const uInt nRows = in.nRow();
818
819// Create output Table and reshape number of polarizations
820
821 Bool clear=True;
822 SDMemTable* pTabOut = new SDMemTable(in, clear);
823 SDHeader header = pTabOut->getSDHeader();
824 header.npol = 1;
825 pTabOut->putSDHeader(header);
826
827// Shape of input and output data
828
829 const IPosition& shapeIn = in.rowAsMaskedArray(0).shape();
830 IPosition shapeOut(shapeIn);
831 shapeOut(asap::PolAxis) = 1; // Average all polarizations
832 if (shapeIn(asap::PolAxis)==1) {
833 throw(AipsError("The input has only one polarisation"));
834 }
835//
836 const uInt nChan = shapeIn(asap::ChanAxis);
837 const IPosition vecShapeOut(4,1,1,1,nChan); // A multi-dim form of a Vector shape
838 IPosition start(4), end(4);
839
840// Output arrays
841
842 Array<Float> outData(shapeOut, 0.0);
843 Array<Bool> outMask(shapeOut, True);
844 const IPosition axes(2, asap::PolAxis, asap::ChanAxis); // pol-channel plane
845//
846 const Bool useMask = (mask.nelements() == shapeIn(asap::ChanAxis));
847
848// Loop over rows
849
850 for (uInt iRow=0; iRow<nRows; iRow++) {
851
852// Get data for this row
853
854 MaskedArray<Float> marr(in.rowAsMaskedArray(iRow));
855 Array<Float>& arr = marr.getRWArray();
856 const Array<Bool>& barr = marr.getMask();
857
858// Make iterators to iterate by pol-channel planes
859
860 ReadOnlyArrayIterator<Float> itDataPlane(arr, axes);
861 ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
862
863// Accumulations
864
865 Float fac = 1.0;
866 Vector<Float> vecSum(nChan,0.0);
867
868// Iterate through data by pol-channel planes
869
870 while (!itDataPlane.pastEnd()) {
871
872// Iterate through plane by polarization and accumulate Vectors
873
874 Vector<Float> t1(nChan); t1 = 0.0;
875 Vector<Bool> t2(nChan); t2 = True;
876 MaskedArray<Float> vecSum(t1,t2);
877 Float norm = 0.0;
878 {
879 ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1);
880 ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
881 while (!itDataVec.pastEnd()) {
882
883// Create MA of data & mask (optionally including OTF mask) and get variance for this spectrum
884
885 if (useMask) {
886 const MaskedArray<Float> spec(itDataVec.vector(),mask&&itMaskVec.vector());
887 if (wtType==VAR) fac = 1.0 / variance(spec);
888 } else {
889 const MaskedArray<Float> spec(itDataVec.vector(),itMaskVec.vector());
890 if (wtType==VAR) fac = 1.0 / variance(spec);
891 }
892
893// Normalize spectrum (without OTF mask) and accumulate
894
895 const MaskedArray<Float> spec(fac*itDataVec.vector(), itMaskVec.vector());
896 vecSum += spec;
897 norm += fac;
898
899// Next
900
901 itDataVec.next();
902 itMaskVec.next();
903 }
904 }
905
906// Normalize summed spectrum
907
908 vecSum /= norm;
909
910// FInd position in input data array. We are iterating by pol-channel
911// plane so all that will change is beam and IF and that's what we want.
912
913 IPosition pos = itDataPlane.pos();
914
915// Write out data. This is a bit messy. We have to reform the Vector
916// accumulator into an Array of shape (1,1,1,nChan)
917
918 start = pos;
919 end = pos;
920 end(asap::ChanAxis) = nChan-1;
921 outData(start,end) = vecSum.getArray().reform(vecShapeOut);
922 outMask(start,end) = vecSum.getMask().reform(vecShapeOut);
923
924// Step to next beam/IF combination
925
926 itDataPlane.next();
927 itMaskPlane.next();
928 }
929
930// Generate output container and write it to output table
931
932 SDContainer sc = in.getSDContainer();
933 sc.resize(shapeOut);
934//
935 putDataInSDC(sc, outData, outMask);
936 pTabOut->putSDContainer(sc);
937 }
938
939// Set polarization cursor to 0
940
941 pTabOut->setPol(0);
942//
943 return pTabOut;
944}
945
946
947SDMemTable* SDMath::smooth(const SDMemTable& in,
948 const casa::String& kernelType,
949 casa::Float width, Bool doAll) const
950//
951// Should smooth TSys as well
952//
953{
954
955// Number of channels
956
957 const uInt nChan = in.nChan();
958
959// Generate Kernel
960
961 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernelType);
962 Vector<Float> kernel = VectorKernel::make(type, width, nChan, True, False);
963
964// Generate Convolver
965
966 IPosition shape(1,nChan);
967 Convolver<Float> conv(kernel, shape);
968
969// New Table
970
971 SDMemTable* pTabOut = new SDMemTable(in,True);
972
973// Output Vectors
974
975 Vector<Float> valuesOut(nChan);
976 Vector<Bool> maskOut(nChan);
977
978// Get data slice bounds
979
980 IPosition start, end;
981 setCursorSlice (start, end, doAll, in);
982
983// Loop over rows in Table
984
985 for (uInt ri=0; ri < in.nRow(); ++ri) {
986
987// Get slice of data
988
989 MaskedArray<Float> dataIn = in.rowAsMaskedArray(ri);
990
991// Deconstruct and get slices which reference these arrays
992
993 Array<Float> valuesIn = dataIn.getArray();
994 Array<Bool> maskIn = dataIn.getMask();
995//
996 Array<Float> valuesIn2 = valuesIn(start,end); // ref to valuesIn
997 Array<Bool> maskIn2 = maskIn(start,end);
998
999// Iterate through by spectra
1000
1001 VectorIterator<Float> itValues(valuesIn2, asap::ChanAxis);
1002 VectorIterator<Bool> itMask(maskIn2, asap::ChanAxis);
1003 while (!itValues.pastEnd()) {
1004
1005// Smooth
1006
1007 if (kernelType==VectorKernel::HANNING) {
1008 mathutil::hanning(valuesOut, maskOut, itValues.vector(), itMask.vector());
1009 itMask.vector() = maskOut;
1010 } else {
1011 mathutil::replaceMaskByZero(itValues.vector(), itMask.vector());
1012 conv.linearConv(valuesOut, itValues.vector());
1013 }
1014//
1015 itValues.vector() = valuesOut;
1016//
1017 itValues.next();
1018 itMask.next();
1019 }
1020
1021// Create and put back
1022
1023 SDContainer sc = in.getSDContainer(ri);
1024 putDataInSDC(sc, valuesIn, maskIn);
1025//
1026 pTabOut->putSDContainer(sc);
1027 }
1028//
1029 return pTabOut;
1030}
1031
1032
1033
1034SDMemTable* SDMath::convertFlux(const SDMemTable& in, Float D, Float etaAp,
1035 Float JyPerK, Bool doAll) const
1036//
1037// etaAp = aperture efficiency (-1 means find)
1038// D = geometric diameter (m) (-1 means find)
1039// JyPerK
1040//
1041{
1042 SDHeader sh = in.getSDHeader();
1043 SDMemTable* pTabOut = new SDMemTable(in, True);
1044
1045// Find out how to convert values into Jy and K (e.g. units might be mJy or mK)
1046// Also automatically find out what we are converting to according to the
1047// flux unit
1048
1049 Unit fluxUnit(sh.fluxunit);
1050 Unit K(String("K"));
1051 Unit JY(String("Jy"));
1052//
1053 Bool toKelvin = True;
1054 Double cFac = 1.0;
1055 if (fluxUnit==JY) {
1056 cout << "Converting to K" << endl;
1057//
1058 Quantum<Double> t(1.0,fluxUnit);
1059 Quantum<Double> t2 = t.get(JY);
1060 cFac = (t2 / t).getValue(); // value to Jy
1061//
1062 toKelvin = True;
1063 sh.fluxunit = "K";
1064 } else if (fluxUnit==K) {
1065 cout << "Converting to Jy" << endl;
1066//
1067 Quantum<Double> t(1.0,fluxUnit);
1068 Quantum<Double> t2 = t.get(K);
1069 cFac = (t2 / t).getValue(); // value to K
1070//
1071 toKelvin = False;
1072 sh.fluxunit = "Jy";
1073 } else {
1074 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
1075 }
1076 pTabOut->putSDHeader(sh);
1077
1078// Make sure input values are converted to either Jy or K first...
1079
1080 Float factor = cFac;
1081
1082// Select method
1083
1084 if (JyPerK>0.0) {
1085 factor *= JyPerK;
1086 if (toKelvin) factor = 1.0 / JyPerK;
1087//
1088 cout << "Jy/K = " << JyPerK << endl;
1089 Vector<Float> factors(in.nRow(), factor);
1090 scaleByVector(pTabOut, in, doAll, factors, False);
1091 } else if (etaAp>0.0) {
1092 Bool throwIt = True;
1093 Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
1094 SDAttr sda;
1095 if (D < 0) D = sda.diameter(inst);
1096 Float JyPerK = SDAttr::findJyPerK (etaAp,D);
1097 cout << "Jy/K = " << JyPerK << endl;
1098 factor *= JyPerK;
1099 if (toKelvin) {
1100 factor = 1.0 / factor;
1101 }
1102//
1103 Vector<Float> factors(in.nRow(), factor);
1104 scaleByVector(pTabOut, in, doAll, factors, False);
1105 } else {
1106
1107// OK now we must deal with automatic look up of values.
1108// We must also deal with the fact that the factors need
1109// to be computed per IF and may be different and may
1110// change per integration.
1111
1112 cout << "Looking up conversion factors" << endl;
1113 convertBrightnessUnits (pTabOut, in, toKelvin, cFac, doAll);
1114 }
1115//
1116 return pTabOut;
1117}
1118
1119
1120
1121
1122
1123SDMemTable* SDMath::gainElevation(const SDMemTable& in,
1124 const Vector<Float>& coeffs,
1125 const String& fileName,
1126 const String& methodStr, Bool doAll) const
1127{
1128
1129// Get header and clone output table
1130
1131 SDHeader sh = in.getSDHeader();
1132 SDMemTable* pTabOut = new SDMemTable(in, True);
1133
1134// Get elevation data from SDMemTable and convert to degrees
1135
1136 const Table& tab = in.table();
1137 ROScalarColumn<Float> elev(tab, "ELEVATION");
1138 Vector<Float> x = elev.getColumn();
1139 x *= Float(180 / C::pi); // Degrees
1140//
1141 const uInt nC = coeffs.nelements();
1142 if (fileName.length()>0 && nC>0) {
1143 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
1144 }
1145
1146// Correct
1147
1148 if (nC>0 || fileName.length()==0) {
1149
1150// Find instrument
1151
1152 Bool throwIt = True;
1153 Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
1154
1155// Set polynomial
1156
1157 Polynomial<Float>* pPoly = 0;
1158 Vector<Float> coeff;
1159 String msg;
1160 if (nC>0) {
1161 pPoly = new Polynomial<Float>(nC);
1162 coeff = coeffs;
1163 msg = String("user");
1164 } else {
1165 SDAttr sdAttr;
1166 coeff = sdAttr.gainElevationPoly(inst);
1167 pPoly = new Polynomial<Float>(3);
1168 msg = String("built in");
1169 }
1170//
1171 if (coeff.nelements()>0) {
1172 pPoly->setCoefficients(coeff);
1173 } else {
1174 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
1175 }
1176//
1177 cout << "Making polynomial correction with " << msg << " coefficients" << endl;
1178 const uInt nRow = in.nRow();
1179 Vector<Float> factor(nRow);
1180 for (uInt i=0; i<nRow; i++) {
1181 factor[i] = 1.0 / (*pPoly)(x[i]);
1182 }
1183 delete pPoly;
1184//
1185 scaleByVector (pTabOut, in, doAll, factor, True);
1186 } else {
1187
1188// Indicate which columns to read from ascii file
1189
1190 String col0("ELEVATION");
1191 String col1("FACTOR");
1192
1193// Read and correct
1194
1195 cout << "Making correction from ascii Table" << endl;
1196 scaleFromAsciiTable (pTabOut, in, fileName, col0, col1,
1197 methodStr, doAll, x, True);
1198 }
1199//
1200 return pTabOut;
1201}
1202
1203
1204
1205SDMemTable* SDMath::opacity(const SDMemTable& in, Float tau, Bool doAll) const
1206{
1207
1208// Get header and clone output table
1209
1210 SDHeader sh = in.getSDHeader();
1211 SDMemTable* pTabOut = new SDMemTable(in, True);
1212
1213// Get elevation data from SDMemTable and convert to degrees
1214
1215 const Table& tab = in.table();
1216 ROScalarColumn<Float> elev(tab, "ELEVATION");
1217 Vector<Float> zDist = elev.getColumn();
1218 zDist = Float(C::pi_2) - zDist;
1219
1220// Generate correction factor
1221
1222 const uInt nRow = in.nRow();
1223 Vector<Float> factor(nRow);
1224 Vector<Float> factor2(nRow);
1225 for (uInt i=0; i<nRow; i++) {
1226 factor[i] = exp(tau)/cos(zDist[i]);
1227 }
1228
1229// Correct
1230
1231 scaleByVector (pTabOut, in, doAll, factor, True);
1232//
1233 return pTabOut;
1234}
1235
1236
1237void SDMath::rotateXYPhase(SDMemTable& in, Float value, Bool doAll)
1238//
1239// phase in degrees
1240// assumes linear correlations
1241//
1242{
1243 if (in.nPol() != 4) {
1244 throw(AipsError("You must have 4 polarizations to run this function"));
1245 }
1246//
1247 SDHeader sh = in.getSDHeader();
1248 Instrument inst = SDAttr::convertInstrument (sh.antennaname, False);
1249 SDAttr sdAtt;
1250 if (sdAtt.feedPolType(inst) != LINEAR) {
1251 throw(AipsError("Only linear polarizations are supported"));
1252 }
1253//
1254 const Table& tabIn = in.table();
1255 ArrayColumn<Float> specCol(tabIn,"SPECTRA");
1256 IPosition start(asap::nAxes,0);
1257 IPosition end(asap::nAxes);
1258
1259// Set cursor slice. Assumes shape the same for all rows
1260
1261 setCursorSlice (start, end, doAll, in);
1262 IPosition start3(start);
1263 start3(asap::PolAxis) = 2; // Real(XY)
1264 IPosition end3(end);
1265 end3(asap::PolAxis) = 2;
1266//
1267 IPosition start4(start);
1268 start4(asap::PolAxis) = 3; // Imag (XY)
1269 IPosition end4(end);
1270 end4(asap::PolAxis) = 3;
1271//
1272 uInt nRow = in.nRow();
1273 Array<Float> data;
1274 for (uInt i=0; i<nRow;++i) {
1275 specCol.get(i,data);
1276 IPosition shape = data.shape();
1277
1278// Get polarization slice references
1279
1280 Array<Float> C3 = data(start3,end3);
1281 Array<Float> C4 = data(start4,end4);
1282
1283// Rotate
1284
1285 SDPolUtil::rotatePhase(C3, C4, value);
1286
1287// Put
1288
1289 specCol.put(i,data);
1290 }
1291}
1292
1293
1294
1295void SDMath::rotateLinPolPhase(SDMemTable& in, Float value, Bool doAll)
1296//
1297// phase in degrees
1298// assumes linear correlations
1299//
1300{
1301 if (in.nPol() != 4) {
1302 throw(AipsError("You must have 4 polarizations to run this function"));
1303 }
1304//
1305 SDHeader sh = in.getSDHeader();
1306 Instrument inst = SDAttr::convertInstrument (sh.antennaname, False);
1307 SDAttr sdAtt;
1308 if (sdAtt.feedPolType(inst) != LINEAR) {
1309 throw(AipsError("Only linear polarizations are supported"));
1310 }
1311//
1312 const Table& tabIn = in.table();
1313 ArrayColumn<Float> specCol(tabIn,"SPECTRA");
1314 ROArrayColumn<Float> stokesCol(tabIn,"STOKES");
1315 IPosition start(asap::nAxes,0);
1316 IPosition end(asap::nAxes);
1317
1318// Set cursor slice. Assumes shape the same for all rows
1319
1320 setCursorSlice (start, end, doAll, in);
1321//
1322 IPosition start1(start);
1323 start1(asap::PolAxis) = 0; // C1 (XX)
1324 IPosition end1(end);
1325 end1(asap::PolAxis) = 0;
1326//
1327 IPosition start2(start);
1328 start2(asap::PolAxis) = 1; // C2 (YY)
1329 IPosition end2(end);
1330 end2(asap::PolAxis) = 1;
1331//
1332 IPosition start3(start);
1333 start3(asap::PolAxis) = 2; // C3 ( Real(XY) )
1334 IPosition end3(end);
1335 end3(asap::PolAxis) = 2;
1336//
1337 IPosition startI(start);
1338 startI(asap::PolAxis) = 0; // I
1339 IPosition endI(end);
1340 endI(asap::PolAxis) = 0;
1341//
1342 IPosition startQ(start);
1343 startQ(asap::PolAxis) = 1; // Q
1344 IPosition endQ(end);
1345 endQ(asap::PolAxis) = 1;
1346//
1347 IPosition startU(start);
1348 startU(asap::PolAxis) = 2; // U
1349 IPosition endU(end);
1350 endU(asap::PolAxis) = 2;
1351
1352//
1353 uInt nRow = in.nRow();
1354 Array<Float> data, stokes;
1355 for (uInt i=0; i<nRow;++i) {
1356 specCol.get(i,data);
1357 stokesCol.get(i,stokes);
1358 IPosition shape = data.shape();
1359
1360// Get linear polarization slice references
1361
1362 Array<Float> C1 = data(start1,end1);
1363 Array<Float> C2 = data(start2,end2);
1364 Array<Float> C3 = data(start3,end3);
1365
1366// Get STokes slice references
1367
1368 Array<Float> I = stokes(startI,endI);
1369 Array<Float> Q = stokes(startQ,endQ);
1370 Array<Float> U = stokes(startU,endU);
1371
1372// Rotate
1373
1374 SDPolUtil::rotateLinPolPhase(C1, C2, C3, I, Q, U, value);
1375
1376// Put
1377
1378 specCol.put(i,data);
1379 }
1380}
1381
1382// 'private' functions
1383
1384void SDMath::convertBrightnessUnits (SDMemTable* pTabOut, const SDMemTable& in,
1385 Bool toKelvin, Float cFac, Bool doAll) const
1386{
1387
1388// Get header
1389
1390 SDHeader sh = in.getSDHeader();
1391 const uInt nChan = sh.nchan;
1392
1393// Get instrument
1394
1395 Bool throwIt = True;
1396 Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
1397
1398// Get Diameter (m)
1399
1400 SDAttr sdAtt;
1401
1402// Get epoch of first row
1403
1404 MEpoch dateObs = in.getEpoch(0);
1405
1406// Generate a Vector of correction factors. One per FreqID
1407
1408 SDFrequencyTable sdft = in.getSDFreqTable();
1409 Vector<uInt> freqIDs;
1410//
1411 Vector<Float> freqs(sdft.length());
1412 for (uInt i=0; i<sdft.length(); i++) {
1413 freqs(i) = (nChan/2 - sdft.referencePixel(i))*sdft.increment(i) + sdft.referenceValue(i);
1414 }
1415//
1416 Vector<Float> JyPerK = sdAtt.JyPerK(inst, dateObs, freqs);
1417 cout << "Jy/K = " << JyPerK << endl;
1418 Vector<Float> factors = cFac * JyPerK;
1419 if (toKelvin) factors = Float(1.0) / factors;
1420
1421// Get data slice bounds
1422
1423 IPosition start, end;
1424 setCursorSlice (start, end, doAll, in);
1425 const uInt ifAxis = in.getIF();
1426
1427// Iteration axes
1428
1429 IPosition axes(asap::nAxes-1,0);
1430 for (uInt i=0,j=0; i<asap::nAxes; i++) {
1431 if (i!=asap::IFAxis) {
1432 axes(j++) = i;
1433 }
1434 }
1435
1436// Loop over rows and apply correction factor
1437
1438 Float factor = 1.0;
1439 const uInt axis = asap::ChanAxis;
1440 for (uInt i=0; i < in.nRow(); ++i) {
1441
1442// Get data
1443
1444 MaskedArray<Float> dataIn = in.rowAsMaskedArray(i);
1445 Array<Float>& values = dataIn.getRWArray(); // Ref to dataIn
1446 Array<Float> values2 = values(start,end); // Ref to values to dataIn
1447
1448// Get SDCOntainer
1449
1450 SDContainer sc = in.getSDContainer(i);
1451
1452// Get FreqIDs
1453
1454 freqIDs = sc.getFreqMap();
1455
1456// Now the conversion factor depends only upon frequency
1457// So we need to iterate through by IF only giving
1458// us BEAM/POL/CHAN cubes
1459
1460 ArrayIterator<Float> itIn(values2, axes);
1461 uInt ax = 0;
1462 while (!itIn.pastEnd()) {
1463 itIn.array() *= factors(freqIDs(ax)); // Writes back to dataIn
1464 itIn.next();
1465 }
1466
1467// Write out
1468
1469 putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
1470//
1471 pTabOut->putSDContainer(sc);
1472 }
1473}
1474
1475
1476
1477SDMemTable* SDMath::frequencyAlign (const SDMemTable& in,
1478 MFrequency::Types freqSystem,
1479 const String& refTime,
1480 const String& methodStr,
1481 Bool perFreqID) const
1482{
1483// Get Header
1484
1485 SDHeader sh = in.getSDHeader();
1486 const uInt nChan = sh.nchan;
1487 const uInt nRows = in.nRow();
1488 const uInt nIF = sh.nif;
1489
1490// Get Table reference
1491
1492 const Table& tabIn = in.table();
1493
1494// Get Columns from Table
1495
1496 ROScalarColumn<Double> mjdCol(tabIn, "TIME");
1497 ROScalarColumn<String> srcCol(tabIn, "SRCNAME");
1498 ROArrayColumn<uInt> fqIDCol(tabIn, "FREQID");
1499 Vector<Double> times = mjdCol.getColumn();
1500
1501// Generate DataDesc table
1502
1503 Matrix<uInt> ddIdx;
1504 SDDataDesc dDesc;
1505 generateDataDescTable (ddIdx, dDesc, nIF, in, tabIn, srcCol, fqIDCol, perFreqID);
1506
1507// Get reference Epoch to time of first row or given String
1508
1509 Unit DAY(String("d"));
1510 MEpoch::Ref epochRef(in.getTimeReference());
1511 MEpoch refEpoch;
1512 if (refTime.length()>0) {
1513 refEpoch = epochFromString(refTime, in.getTimeReference());
1514 } else {
1515 refEpoch = in.getEpoch(0);
1516 }
1517 cout << "Aligning at reference Epoch " << formatEpoch(refEpoch)
1518 << " in frame " << MFrequency::showType(freqSystem) << endl;
1519
1520// Get Reference Position
1521
1522 MPosition refPos = in.getAntennaPosition();
1523
1524// Create FrequencyAligner Block. One FA for each possible
1525// source/freqID (perFreqID=True) or source/IF (perFreqID=False) combination
1526
1527 PtrBlock<FrequencyAligner<Float>* > a(dDesc.length());
1528 generateFrequencyAligners (a, dDesc, in, nChan, freqSystem, refPos,
1529 refEpoch, perFreqID);
1530
1531// Generate and fill output Frequency Table. WHen perFreqID=True, there is one output FreqID
1532// for each entry in the SDDataDesc table. However, in perFreqID=False mode, there may be
1533// some degeneracy, so we need a little translation map
1534
1535 SDFrequencyTable freqTabOut = in.getSDFreqTable();
1536 freqTabOut.setLength(0);
1537 Vector<String> units(1);
1538 units = String("Hz");
1539 Bool linear=True;
1540//
1541 Vector<uInt> ddFQTrans(dDesc.length(),0);
1542 for (uInt i=0; i<dDesc.length(); i++) {
1543
1544// Get Aligned SC in Hz
1545
1546 SpectralCoordinate sC = a[i]->alignedSpectralCoordinate(linear);
1547 sC.setWorldAxisUnits(units);
1548
1549// Add FreqID
1550
1551 uInt idx = freqTabOut.addFrequency(sC.referencePixel()[0],
1552 sC.referenceValue()[0],
1553 sC.increment()[0]);
1554 ddFQTrans(i) = idx; // output FreqID = ddFQTrans(ddIdx)
1555 }
1556
1557// Interpolation method
1558
1559 InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1560 convertInterpString(interp, methodStr);
1561
1562// New output Table
1563
1564 cout << "Create output table" << endl;
1565 SDMemTable* pTabOut = new SDMemTable(in,True);
1566 pTabOut->putSDFreqTable(freqTabOut);
1567
1568// Loop over rows in Table
1569
1570 Bool extrapolate=False;
1571 const IPosition polChanAxes(2, asap::PolAxis, asap::ChanAxis);
1572 Bool useCachedAbcissa = False;
1573 Bool first = True;
1574 Bool ok;
1575 Vector<Float> yOut;
1576 Vector<Bool> maskOut;
1577 Vector<uInt> freqID(nIF);
1578 uInt ifIdx, faIdx;
1579 Vector<Double> xIn;
1580//
1581 for (uInt iRow=0; iRow<nRows; ++iRow) {
1582 if (iRow%10==0) {
1583 cout << "Processing row " << iRow << endl;
1584 }
1585
1586// Get EPoch
1587
1588 Quantum<Double> tQ2(times[iRow],DAY);
1589 MVEpoch mv2(tQ2);
1590 MEpoch epoch(mv2, epochRef);
1591
1592// Get copy of data
1593
1594 const MaskedArray<Float>& mArrIn(in.rowAsMaskedArray(iRow));
1595 Array<Float> values = mArrIn.getArray();
1596 Array<Bool> mask = mArrIn.getMask();
1597
1598// For each row, the Frequency abcissa will be the same regardless
1599// of polarization. For all other axes (IF and BEAM) the abcissa
1600// will change. So we iterate through the data by pol-chan planes
1601// to mimimize the work. Probably won't work for multiple beams
1602// at this point.
1603
1604 ArrayIterator<Float> itValuesPlane(values, polChanAxes);
1605 ArrayIterator<Bool> itMaskPlane(mask, polChanAxes);
1606 while (!itValuesPlane.pastEnd()) {
1607
1608// Find the IF index and then the FA PtrBlock index
1609
1610 const IPosition& pos = itValuesPlane.pos();
1611 ifIdx = pos(asap::IFAxis);
1612 faIdx = ddIdx(iRow,ifIdx);
1613
1614// Generate abcissa for perIF. Could cache this in a Matrix
1615// on a per scan basis. Pretty expensive doing it for every row.
1616
1617 if (!perFreqID) {
1618 xIn.resize(nChan);
1619 uInt fqID = dDesc.secID(ddIdx(iRow,ifIdx));
1620 SpectralCoordinate sC = in.getSpectralCoordinate(fqID);
1621 Double w;
1622 for (uInt i=0; i<nChan; i++) {
1623 sC.toWorld(w,Double(i));
1624 xIn[i] = w;
1625 }
1626 }
1627//
1628 VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1);
1629 VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
1630
1631// Iterate through the plane by vector and align
1632
1633 first = True;
1634 useCachedAbcissa=False;
1635 while (!itValuesVec.pastEnd()) {
1636 if (perFreqID) {
1637 ok = a[faIdx]->align (yOut, maskOut, itValuesVec.vector(),
1638 itMaskVec.vector(), epoch, useCachedAbcissa,
1639 interp, extrapolate);
1640 } else {
1641 ok = a[faIdx]->align (yOut, maskOut, xIn, itValuesVec.vector(),
1642 itMaskVec.vector(), epoch, useCachedAbcissa,
1643 interp, extrapolate);
1644 }
1645//
1646 itValuesVec.vector() = yOut;
1647 itMaskVec.vector() = maskOut;
1648//
1649 itValuesVec.next();
1650 itMaskVec.next();
1651//
1652 if (first) {
1653 useCachedAbcissa = True;
1654 first = False;
1655 }
1656 }
1657//
1658 itValuesPlane.next();
1659 itMaskPlane.next();
1660 }
1661
1662// Create SDContainer and put back
1663
1664 SDContainer sc = in.getSDContainer(iRow);
1665 putDataInSDC(sc, values, mask);
1666
1667// Set output FreqIDs
1668
1669 for (uInt i=0; i<nIF; i++) {
1670 uInt idx = ddIdx(iRow,i); // Index into SDDataDesc table
1671 freqID(i) = ddFQTrans(idx); // FreqID in output FQ table
1672 }
1673 sc.putFreqMap(freqID);
1674//
1675 pTabOut->putSDContainer(sc);
1676 }
1677
1678// Now we must set the base and extra frames to the
1679// input frame
1680
1681 std::vector<string> info = pTabOut->getCoordInfo();
1682 info[1] = MFrequency::showType(freqSystem); // Conversion frame
1683 info[3] = info[1]; // Base frame
1684 pTabOut->setCoordInfo(info);
1685
1686// Clean up PointerBlock
1687
1688 for (uInt i=0; i<a.nelements(); i++) delete a[i];
1689//
1690 return pTabOut;
1691}
1692
1693
1694void SDMath::fillSDC(SDContainer& sc,
1695 const Array<Bool>& mask,
1696 const Array<Float>& data,
1697 const Array<Float>& tSys,
1698 Int scanID, Double timeStamp,
1699 Double interval, const String& sourceName,
1700 const Vector<uInt>& freqID) const
1701{
1702// Data and mask
1703
1704 putDataInSDC(sc, data, mask);
1705
1706// TSys
1707
1708 sc.putTsys(tSys);
1709
1710// Time things
1711
1712 sc.timestamp = timeStamp;
1713 sc.interval = interval;
1714 sc.scanid = scanID;
1715//
1716 sc.sourcename = sourceName;
1717 sc.putFreqMap(freqID);
1718}
1719
1720void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum,
1721 MaskedArray<Float>& sum, Array<Float>& sumSq,
1722 Array<Float>& nPts, Array<Float>& tSysSum,
1723 Array<Float>& tSysSqSum,
1724 const Array<Float>& tSys, const Array<Float>& nInc,
1725 const Vector<Bool>& mask, Double time, Double interval,
1726 const Block<CountedPtr<SDMemTable> >& in,
1727 uInt iTab, uInt iRow, uInt axis,
1728 uInt nAxesSub, Bool useMask,
1729 WeightType wtType) const
1730{
1731
1732// Get data
1733
1734 MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
1735 Array<Float>& valuesIn = dataIn.getRWArray(); // writable reference
1736 const Array<Bool>& maskIn = dataIn.getMask(); // RO reference
1737//
1738 if (wtType==NONE) {
1739 const MaskedArray<Float> n(nInc,dataIn.getMask());
1740 nPts += n; // Only accumulates where mask==T
1741 } else if (wtType==TINT) {
1742
1743// We are weighting the data by integration time.
1744
1745 valuesIn *= Float(interval);
1746
1747 } else if (wtType==VAR) {
1748
1749// We are going to average the data, weighted by the noise for each pol, beam and IF.
1750// So therefore we need to iterate through by spectrum (axis 3)
1751
1752 VectorIterator<Float> itData(valuesIn, axis);
1753 ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
1754 Float fac = 1.0;
1755 IPosition pos(nAxesSub,0);
1756//
1757 while (!itData.pastEnd()) {
1758
1759// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
1760
1761 if (useMask) {
1762 MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
1763 fac = 1.0/variance(tmp);
1764 } else {
1765 MaskedArray<Float> tmp(itData.vector(),itMask.vector());
1766 fac = 1.0/variance(tmp);
1767 }
1768
1769// Scale data
1770
1771 itData.vector() *= fac; // Writes back into 'dataIn'
1772//
1773// Accumulate variance per if/pol/beam averaged over spectrum
1774// This method to get pos2 from itData.pos() is only valid
1775// because the spectral axis is the last one (so we can just
1776// copy the first nAXesSub positions out)
1777
1778 pos = itData.pos().getFirst(nAxesSub);
1779 sumSq(pos) += fac;
1780//
1781 itData.next();
1782 itMask.next();
1783 }
1784 } else if (wtType==TSYS) {
1785
1786// We are going to average the data, weighted by 1/Tsys**2 for each pol, beam and IF.
1787// So therefore we need to iterate through by spectrum (axis 3). Although
1788// Tsys is stored as a vector of length nChan, the values are replicated.
1789// We will take a short cut and just use the value from the first channel
1790// for now.
1791//
1792 VectorIterator<Float> itData(valuesIn, axis);
1793 ReadOnlyVectorIterator<Float> itTSys(tSys, axis);
1794 Float fac = 1.0;
1795 IPosition pos(nAxesSub,0);
1796//
1797 while (!itData.pastEnd()) {
1798 Float t = itTSys.vector()[0];
1799 fac = 1.0/t/t;
1800
1801// Scale data
1802
1803 itData.vector() *= fac; // Writes back into 'dataIn'
1804//
1805// Accumulate Tsys per if/pol/beam averaged over spectrum
1806// This method to get pos2 from itData.pos() is only valid
1807// because the spectral axis is the last one (so we can just
1808// copy the first nAXesSub positions out)
1809
1810 pos = itData.pos().getFirst(nAxesSub);
1811 tSysSqSum(pos) += fac;
1812//
1813 itData.next();
1814 itTSys.next();
1815 }
1816 }
1817
1818// Accumulate sum of (possibly scaled) data
1819
1820 sum += dataIn;
1821
1822// Accumulate Tsys, time, and interval
1823
1824 tSysSum += tSys;
1825 timeSum += time;
1826 intSum += interval;
1827 nAccum += 1;
1828}
1829
1830
1831void SDMath::normalize(MaskedArray<Float>& sum,
1832 const Array<Float>& sumSq,
1833 const Array<Float>& tSysSqSum,
1834 const Array<Float>& nPts,
1835 Double intSum,
1836 WeightType wtType, Int axis,
1837 Int nAxesSub) const
1838{
1839 IPosition pos2(nAxesSub,0);
1840//
1841 if (wtType==NONE) {
1842
1843// We just average by the number of points accumulated.
1844// We need to make a MA out of nPts so that no divide by
1845// zeros occur
1846
1847 MaskedArray<Float> t(nPts, (nPts>Float(0.0)));
1848 sum /= t;
1849 } else if (wtType==TINT) {
1850
1851// Average by sum of Tint
1852
1853 sum /= Float(intSum);
1854 } else if (wtType==VAR) {
1855
1856// Normalize each spectrum by sum(1/var) where the variance
1857// is worked out for each spectrum
1858
1859 Array<Float>& data = sum.getRWArray();
1860 VectorIterator<Float> itData(data, axis);
1861 while (!itData.pastEnd()) {
1862 pos2 = itData.pos().getFirst(nAxesSub);
1863 itData.vector() /= sumSq(pos2);
1864 itData.next();
1865 }
1866 } else if (wtType==TSYS) {
1867
1868// Normalize each spectrum by sum(1/Tsys**2) where the pseudo
1869// replication over channel for Tsys has been dropped.
1870
1871 Array<Float>& data = sum.getRWArray();
1872 VectorIterator<Float> itData(data, axis);
1873 while (!itData.pastEnd()) {
1874 pos2 = itData.pos().getFirst(nAxesSub);
1875 itData.vector() /= tSysSqSum(pos2);
1876 itData.next();
1877 }
1878 }
1879}
1880
1881
1882
1883
1884void SDMath::setCursorSlice (IPosition& start, IPosition& end, Bool doAll, const SDMemTable& in) const
1885{
1886 const uInt nDim = asap::nAxes;
1887 DebugAssert(nDim==4,AipsError);
1888//
1889 start.resize(nDim);
1890 end.resize(nDim);
1891 if (doAll) {
1892 start = 0;
1893 end(0) = in.nBeam()-1;
1894 end(1) = in.nIF()-1;
1895 end(2) = in.nPol()-1;
1896 end(3) = in.nChan()-1;
1897 } else {
1898 start(0) = in.getBeam();
1899 end(0) = start(0);
1900//
1901 start(1) = in.getIF();
1902 end(1) = start(1);
1903//
1904 start(2) = in.getPol();
1905 end(2) = start(2);
1906//
1907 start(3) = 0;
1908 end(3) = in.nChan()-1;
1909 }
1910}
1911
1912
1913void SDMath::convertWeightString(WeightType& wtType, const String& weightStr,
1914 Bool listType) const
1915{
1916 String tStr(weightStr);
1917 tStr.upcase();
1918 String msg;
1919 if (tStr.contains(String("NONE"))) {
1920 wtType = NONE;
1921 msg = String("Weighting type selected : None");
1922 } else if (tStr.contains(String("VAR"))) {
1923 wtType = VAR;
1924 msg = String("Weighting type selected : Variance");
1925 } else if (tStr.contains(String("TINT"))) {
1926 wtType = TINT;
1927 msg = String("Weighting type selected : Tint");
1928 } else if (tStr.contains(String("TSYS"))) {
1929 wtType = TSYS;
1930 msg = String("Weighting type selected : Tsys");
1931 } else {
1932 msg = String("Weighting type selected : None");
1933 throw(AipsError("Unrecognized weighting type"));
1934 }
1935//
1936 if (listType) cout << msg << endl;
1937}
1938
1939
1940void SDMath::convertInterpString(casa::InterpolateArray1D<Double,Float>::InterpolationMethod& type,
1941 const casa::String& interp) const
1942{
1943 String tStr(interp);
1944 tStr.upcase();
1945 if (tStr.contains(String("NEAR"))) {
1946 type = InterpolateArray1D<Double,Float>::nearestNeighbour;
1947 } else if (tStr.contains(String("LIN"))) {
1948 type = InterpolateArray1D<Double,Float>::linear;
1949 } else if (tStr.contains(String("CUB"))) {
1950 type = InterpolateArray1D<Double,Float>::cubic;
1951 } else if (tStr.contains(String("SPL"))) {
1952 type = InterpolateArray1D<Double,Float>::spline;
1953 } else {
1954 throw(AipsError("Unrecognized interpolation type"));
1955 }
1956}
1957
1958void SDMath::putDataInSDC(SDContainer& sc, const Array<Float>& data,
1959 const Array<Bool>& mask) const
1960{
1961 sc.putSpectrum(data);
1962//
1963 Array<uChar> outflags(data.shape());
1964 convertArray(outflags,!mask);
1965 sc.putFlags(outflags);
1966}
1967
1968Table SDMath::readAsciiFile (const String& fileName) const
1969{
1970 String formatString;
1971 Table tbl = readAsciiTable (formatString, Table::Memory, fileName, "", "", False);
1972 return tbl;
1973}
1974
1975
1976
1977void SDMath::scaleFromAsciiTable(SDMemTable* pTabOut,
1978 const SDMemTable& in, const String& fileName,
1979 const String& col0, const String& col1,
1980 const String& methodStr, Bool doAll,
1981 const Vector<Float>& xOut, Bool doTSys) const
1982{
1983
1984// Read gain-elevation ascii file data into a Table.
1985
1986 Table geTable = readAsciiFile (fileName);
1987//
1988 scaleFromTable (pTabOut, in, geTable, col0, col1, methodStr, doAll, xOut, doTSys);
1989}
1990
1991void SDMath::scaleFromTable(SDMemTable* pTabOut, const SDMemTable& in,
1992 const Table& tTable, const String& col0,
1993 const String& col1,
1994 const String& methodStr, Bool doAll,
1995 const Vector<Float>& xOut, Bool doTsys) const
1996{
1997
1998// Get data from Table
1999
2000 ROScalarColumn<Float> geElCol(tTable, col0);
2001 ROScalarColumn<Float> geFacCol(tTable, col1);
2002 Vector<Float> xIn = geElCol.getColumn();
2003 Vector<Float> yIn = geFacCol.getColumn();
2004 Vector<Bool> maskIn(xIn.nelements(),True);
2005
2006// Interpolate (and extrapolate) with desired method
2007
2008 InterpolateArray1D<Double,Float>::InterpolationMethod method;
2009 convertInterpString(method, methodStr);
2010 Int intMethod(method);
2011//
2012 Vector<Float> yOut;
2013 Vector<Bool> maskOut;
2014 InterpolateArray1D<Float,Float>::interpolate(yOut, maskOut, xOut,
2015 xIn, yIn, maskIn, intMethod,
2016 True, True);
2017// Apply
2018
2019 scaleByVector(pTabOut, in, doAll, Float(1.0)/yOut, doTsys);
2020}
2021
2022
2023void SDMath::scaleByVector(SDMemTable* pTabOut, const SDMemTable& in,
2024 Bool doAll, const Vector<Float>& factor,
2025 Bool doTSys) const
2026{
2027
2028// Set up data slice
2029
2030 IPosition start, end;
2031 setCursorSlice (start, end, doAll, in);
2032
2033// Get Tsys column
2034
2035 const Table& tIn = in.table();
2036 ArrayColumn<Float> tSysCol(tIn, "TSYS");
2037 Array<Float> tSys;
2038
2039// Loop over rows and apply correction factor
2040
2041 const uInt axis = asap::ChanAxis;
2042 for (uInt i=0; i < in.nRow(); ++i) {
2043
2044// Get data
2045
2046 MaskedArray<Float> dataIn(in.rowAsMaskedArray(i));
2047 MaskedArray<Float> dataIn2 = dataIn(start,end); // reference to dataIn
2048//
2049 if (doTSys) {
2050 tSysCol.get(i, tSys);
2051 Array<Float> tSys2 = tSys(start,end) * factor[i];
2052 tSysCol.put(i, tSys);
2053 }
2054
2055// Apply factor
2056
2057 dataIn2 *= factor[i];
2058
2059// Write out
2060
2061 SDContainer sc = in.getSDContainer(i);
2062 putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
2063//
2064 pTabOut->putSDContainer(sc);
2065 }
2066}
2067
2068
2069
2070
2071void SDMath::generateDataDescTable (Matrix<uInt>& ddIdx,
2072 SDDataDesc& dDesc,
2073 uInt nIF,
2074 const SDMemTable& in,
2075 const Table& tabIn,
2076 const ROScalarColumn<String>& srcCol,
2077 const ROArrayColumn<uInt>& fqIDCol,
2078 Bool perFreqID) const
2079{
2080 const uInt nRows = tabIn.nrow();
2081 ddIdx.resize(nRows,nIF);
2082//
2083 String srcName;
2084 Vector<uInt> freqIDs;
2085 for (uInt iRow=0; iRow<nRows; iRow++) {
2086 srcCol.get(iRow, srcName);
2087 fqIDCol.get(iRow, freqIDs);
2088 const MDirection& dir = in.getDirection(iRow);
2089//
2090 if (perFreqID) {
2091
2092// One entry per source/freqID pair
2093
2094 for (uInt iIF=0; iIF<nIF; iIF++) {
2095 ddIdx(iRow,iIF) = dDesc.addEntry(srcName, freqIDs[iIF], dir, 0);
2096 }
2097 } else {
2098
2099// One entry per source/IF pair. Hang onto the FreqID as well
2100
2101 for (uInt iIF=0; iIF<nIF; iIF++) {
2102 ddIdx(iRow,iIF) = dDesc.addEntry(srcName, iIF, dir, freqIDs[iIF]);
2103 }
2104 }
2105 }
2106}
2107
2108
2109
2110
2111
2112MEpoch SDMath::epochFromString (const String& str, MEpoch::Types timeRef) const
2113{
2114 Quantum<Double> qt;
2115 if (MVTime::read(qt,str)) {
2116 MVEpoch mv(qt);
2117 MEpoch me(mv, timeRef);
2118 return me;
2119 } else {
2120 throw(AipsError("Invalid format for Epoch string"));
2121 }
2122}
2123
2124
2125String SDMath::formatEpoch(const MEpoch& epoch) const
2126{
2127 MVTime mvt(epoch.getValue());
2128 return mvt.string(MVTime::YMD) + String(" (") + epoch.getRefString() + String(")");
2129}
2130
2131
2132
2133void SDMath::generateFrequencyAligners (PtrBlock<FrequencyAligner<Float>* >& a,
2134 const SDDataDesc& dDesc,
2135 const SDMemTable& in, uInt nChan,
2136 MFrequency::Types system,
2137 const MPosition& refPos,
2138 const MEpoch& refEpoch,
2139 Bool perFreqID) const
2140{
2141 for (uInt i=0; i<dDesc.length(); i++) {
2142 uInt ID = dDesc.ID(i);
2143 uInt secID = dDesc.secID(i);
2144 const MDirection& refDir = dDesc.secDir(i);
2145//
2146 if (perFreqID) {
2147
2148// One aligner per source/FreqID pair.
2149
2150 SpectralCoordinate sC = in.getSpectralCoordinate(ID);
2151 a[i] = new FrequencyAligner<Float>(sC, nChan, refEpoch, refDir, refPos, system);
2152 } else {
2153
2154// One aligner per source/IF pair. But we still need the FreqID to
2155// get the right SC. Hence the messing about with the secondary ID
2156
2157 SpectralCoordinate sC = in.getSpectralCoordinate(secID);
2158 a[i] = new FrequencyAligner<Float>(sC, nChan, refEpoch, refDir, refPos, system);
2159 }
2160 }
2161}
2162
2163Vector<uInt> SDMath::getRowRange (const SDMemTable& in) const
2164{
2165 Vector<uInt> range(2);
2166 range[0] = 0;
2167 range[1] = in.nRow()-1;
2168 return range;
2169}
2170
2171
2172Bool SDMath::rowInRange (uInt i, const Vector<uInt>& range) const
2173{
2174 return (i>=range[0] && i<=range[1]);
2175}
Note: See TracBrowser for help on using the repository browser.