source: trunk/src/SDMath.cc@ 533

Last change on this file since 533 was 532, checked in by kil064, 20 years ago

add Tsys weighting to average_pol
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 60.3 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, True);
816
817// Create output Table and reshape number of polarizations
818
819 Bool clear=True;
820 SDMemTable* pTabOut = new SDMemTable(in, clear);
821 SDHeader header = pTabOut->getSDHeader();
822 header.npol = 1;
823 pTabOut->putSDHeader(header);
824//
825 const Table& tabIn = in.table();
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 nRows = in.nRow();
837 const uInt nChan = shapeIn(asap::ChanAxis);
838 AlwaysAssert(asap::nAxes==4,AipsError);
839 const IPosition vecShapeOut(4,1,1,1,nChan); // A multi-dim form of a Vector shape
840 IPosition start(4), end(4);
841
842// Output arrays
843
844 Array<Float> outData(shapeOut, 0.0);
845 Array<Bool> outMask(shapeOut, True);
846 const IPosition axes(2, asap::PolAxis, asap::ChanAxis); // pol-channel plane
847
848// Attach Tsys column if needed
849
850 ROArrayColumn<Float> tSysCol;
851 Array<Float> tSys;
852 if (wtType==TSYS) {
853 tSysCol.attach(tabIn,"TSYS");
854 }
855//
856 const Bool useMask = (mask.nelements() == shapeIn(asap::ChanAxis));
857
858// Loop over rows
859
860 for (uInt iRow=0; iRow<nRows; iRow++) {
861
862// Get data for this row
863
864 MaskedArray<Float> marr(in.rowAsMaskedArray(iRow));
865 Array<Float>& arr = marr.getRWArray();
866 const Array<Bool>& barr = marr.getMask();
867
868// Get Tsys
869
870 if (wtType==TSYS) {
871 tSysCol.get(iRow,tSys);
872 }
873
874// Make iterators to iterate by pol-channel planes
875// The tSys array is empty unless wtType=TSYS so only
876// access the iterator is that is the case
877
878 ReadOnlyArrayIterator<Float> itDataPlane(arr, axes);
879 ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
880 ReadOnlyArrayIterator<Float>* pItTsysPlane = 0;
881 if (wtType==TSYS) pItTsysPlane = new ReadOnlyArrayIterator<Float>(tSys, axes);
882
883// Accumulations
884
885 Float fac = 1.0;
886 Vector<Float> vecSum(nChan,0.0);
887
888// Iterate through data by pol-channel planes
889
890 while (!itDataPlane.pastEnd()) {
891
892// Iterate through plane by polarization and accumulate Vectors
893
894 Vector<Float> t1(nChan); t1 = 0.0;
895 Vector<Bool> t2(nChan); t2 = True;
896 Float tSys = 0.0;
897 MaskedArray<Float> vecSum(t1,t2);
898 Float norm = 0.0;
899 {
900 ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1);
901 ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
902//
903 ReadOnlyVectorIterator<Float>* pItTsysVec = 0;
904 if (wtType==TSYS) {
905 pItTsysVec = new ReadOnlyVectorIterator<Float>(pItTsysPlane->array(), 1);
906 }
907//
908 while (!itDataVec.pastEnd()) {
909
910// Create MA of data & mask (optionally including OTF mask) and get variance for this spectrum
911
912 if (useMask) {
913 const MaskedArray<Float> spec(itDataVec.vector(),mask&&itMaskVec.vector());
914 if (wtType==VAR) {
915 fac = 1.0 / variance(spec);
916 } else if (wtType==TSYS) {
917 tSys = pItTsysVec->vector()[0]; // Drop pseudo channel dependency
918 fac = 1.0 / tSys / tSys;
919 }
920 } else {
921 const MaskedArray<Float> spec(itDataVec.vector(),itMaskVec.vector());
922 if (wtType==VAR) {
923 fac = 1.0 / variance(spec);
924 } else if (wtType==TSYS) {
925 tSys = pItTsysVec->vector()[0]; // Drop pseudo channel dependency
926 fac = 1.0 / tSys / tSys;
927 }
928 }
929
930// Normalize spectrum (without OTF mask) and accumulate
931
932 const MaskedArray<Float> spec(fac*itDataVec.vector(), itMaskVec.vector());
933 vecSum += spec;
934 norm += fac;
935
936// Next
937
938 itDataVec.next();
939 itMaskVec.next();
940 if (wtType==TSYS) pItTsysVec->next();
941 }
942
943// Clean up
944
945 if (pItTsysVec) {
946 delete pItTsysVec;
947 pItTsysVec = 0;
948 }
949 }
950
951// Normalize summed spectrum
952
953 vecSum /= norm;
954
955// FInd position in input data array. We are iterating by pol-channel
956// plane so all that will change is beam and IF and that's what we want.
957
958 IPosition pos = itDataPlane.pos();
959
960// Write out data. This is a bit messy. We have to reform the Vector
961// accumulator into an Array of shape (1,1,1,nChan)
962
963 start = pos;
964 end = pos;
965 end(asap::ChanAxis) = nChan-1;
966 outData(start,end) = vecSum.getArray().reform(vecShapeOut);
967 outMask(start,end) = vecSum.getMask().reform(vecShapeOut);
968
969// Step to next beam/IF combination
970
971 itDataPlane.next();
972 itMaskPlane.next();
973 if (wtType==TSYS) pItTsysPlane->next();
974 }
975
976// Generate output container and write it to output table
977
978 SDContainer sc = in.getSDContainer();
979 sc.resize(shapeOut);
980//
981 putDataInSDC(sc, outData, outMask);
982 pTabOut->putSDContainer(sc);
983//
984 if (wtType==TSYS) {
985 delete pItTsysPlane;
986 pItTsysPlane = 0;
987 }
988 }
989
990// Set polarization cursor to 0
991
992 pTabOut->setPol(0);
993//
994 return pTabOut;
995}
996
997
998SDMemTable* SDMath::smooth(const SDMemTable& in,
999 const casa::String& kernelType,
1000 casa::Float width, Bool doAll) const
1001//
1002// Should smooth TSys as well
1003//
1004{
1005
1006// Number of channels
1007
1008 const uInt nChan = in.nChan();
1009
1010// Generate Kernel
1011
1012 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernelType);
1013 Vector<Float> kernel = VectorKernel::make(type, width, nChan, True, False);
1014
1015// Generate Convolver
1016
1017 IPosition shape(1,nChan);
1018 Convolver<Float> conv(kernel, shape);
1019
1020// New Table
1021
1022 SDMemTable* pTabOut = new SDMemTable(in,True);
1023
1024// Output Vectors
1025
1026 Vector<Float> valuesOut(nChan);
1027 Vector<Bool> maskOut(nChan);
1028
1029// Get data slice bounds
1030
1031 IPosition start, end;
1032 setCursorSlice (start, end, doAll, in);
1033
1034// Loop over rows in Table
1035
1036 for (uInt ri=0; ri < in.nRow(); ++ri) {
1037
1038// Get slice of data
1039
1040 MaskedArray<Float> dataIn = in.rowAsMaskedArray(ri);
1041
1042// Deconstruct and get slices which reference these arrays
1043
1044 Array<Float> valuesIn = dataIn.getArray();
1045 Array<Bool> maskIn = dataIn.getMask();
1046//
1047 Array<Float> valuesIn2 = valuesIn(start,end); // ref to valuesIn
1048 Array<Bool> maskIn2 = maskIn(start,end);
1049
1050// Iterate through by spectra
1051
1052 VectorIterator<Float> itValues(valuesIn2, asap::ChanAxis);
1053 VectorIterator<Bool> itMask(maskIn2, asap::ChanAxis);
1054 while (!itValues.pastEnd()) {
1055
1056// Smooth
1057
1058 if (kernelType==VectorKernel::HANNING) {
1059 mathutil::hanning(valuesOut, maskOut, itValues.vector(), itMask.vector());
1060 itMask.vector() = maskOut;
1061 } else {
1062 mathutil::replaceMaskByZero(itValues.vector(), itMask.vector());
1063 conv.linearConv(valuesOut, itValues.vector());
1064 }
1065//
1066 itValues.vector() = valuesOut;
1067//
1068 itValues.next();
1069 itMask.next();
1070 }
1071
1072// Create and put back
1073
1074 SDContainer sc = in.getSDContainer(ri);
1075 putDataInSDC(sc, valuesIn, maskIn);
1076//
1077 pTabOut->putSDContainer(sc);
1078 }
1079//
1080 return pTabOut;
1081}
1082
1083
1084
1085SDMemTable* SDMath::convertFlux(const SDMemTable& in, Float D, Float etaAp,
1086 Float JyPerK, Bool doAll) const
1087//
1088// etaAp = aperture efficiency (-1 means find)
1089// D = geometric diameter (m) (-1 means find)
1090// JyPerK
1091//
1092{
1093 SDHeader sh = in.getSDHeader();
1094 SDMemTable* pTabOut = new SDMemTable(in, True);
1095
1096// Find out how to convert values into Jy and K (e.g. units might be mJy or mK)
1097// Also automatically find out what we are converting to according to the
1098// flux unit
1099
1100 Unit fluxUnit(sh.fluxunit);
1101 Unit K(String("K"));
1102 Unit JY(String("Jy"));
1103//
1104 Bool toKelvin = True;
1105 Double cFac = 1.0;
1106 if (fluxUnit==JY) {
1107 cout << "Converting to K" << endl;
1108//
1109 Quantum<Double> t(1.0,fluxUnit);
1110 Quantum<Double> t2 = t.get(JY);
1111 cFac = (t2 / t).getValue(); // value to Jy
1112//
1113 toKelvin = True;
1114 sh.fluxunit = "K";
1115 } else if (fluxUnit==K) {
1116 cout << "Converting to Jy" << endl;
1117//
1118 Quantum<Double> t(1.0,fluxUnit);
1119 Quantum<Double> t2 = t.get(K);
1120 cFac = (t2 / t).getValue(); // value to K
1121//
1122 toKelvin = False;
1123 sh.fluxunit = "Jy";
1124 } else {
1125 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
1126 }
1127 pTabOut->putSDHeader(sh);
1128
1129// Make sure input values are converted to either Jy or K first...
1130
1131 Float factor = cFac;
1132
1133// Select method
1134
1135 if (JyPerK>0.0) {
1136 factor *= JyPerK;
1137 if (toKelvin) factor = 1.0 / JyPerK;
1138//
1139 cout << "Jy/K = " << JyPerK << endl;
1140 Vector<Float> factors(in.nRow(), factor);
1141 scaleByVector(pTabOut, in, doAll, factors, False);
1142 } else if (etaAp>0.0) {
1143 Bool throwIt = True;
1144 Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
1145 SDAttr sda;
1146 if (D < 0) D = sda.diameter(inst);
1147 Float JyPerK = SDAttr::findJyPerK (etaAp,D);
1148 cout << "Jy/K = " << JyPerK << endl;
1149 factor *= JyPerK;
1150 if (toKelvin) {
1151 factor = 1.0 / factor;
1152 }
1153//
1154 Vector<Float> factors(in.nRow(), factor);
1155 scaleByVector(pTabOut, in, doAll, factors, False);
1156 } else {
1157
1158// OK now we must deal with automatic look up of values.
1159// We must also deal with the fact that the factors need
1160// to be computed per IF and may be different and may
1161// change per integration.
1162
1163 cout << "Looking up conversion factors" << endl;
1164 convertBrightnessUnits (pTabOut, in, toKelvin, cFac, doAll);
1165 }
1166//
1167 return pTabOut;
1168}
1169
1170
1171
1172
1173
1174SDMemTable* SDMath::gainElevation(const SDMemTable& in,
1175 const Vector<Float>& coeffs,
1176 const String& fileName,
1177 const String& methodStr, Bool doAll) const
1178{
1179
1180// Get header and clone output table
1181
1182 SDHeader sh = in.getSDHeader();
1183 SDMemTable* pTabOut = new SDMemTable(in, True);
1184
1185// Get elevation data from SDMemTable and convert to degrees
1186
1187 const Table& tab = in.table();
1188 ROScalarColumn<Float> elev(tab, "ELEVATION");
1189 Vector<Float> x = elev.getColumn();
1190 x *= Float(180 / C::pi); // Degrees
1191//
1192 const uInt nC = coeffs.nelements();
1193 if (fileName.length()>0 && nC>0) {
1194 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
1195 }
1196
1197// Correct
1198
1199 if (nC>0 || fileName.length()==0) {
1200
1201// Find instrument
1202
1203 Bool throwIt = True;
1204 Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
1205
1206// Set polynomial
1207
1208 Polynomial<Float>* pPoly = 0;
1209 Vector<Float> coeff;
1210 String msg;
1211 if (nC>0) {
1212 pPoly = new Polynomial<Float>(nC);
1213 coeff = coeffs;
1214 msg = String("user");
1215 } else {
1216 SDAttr sdAttr;
1217 coeff = sdAttr.gainElevationPoly(inst);
1218 pPoly = new Polynomial<Float>(3);
1219 msg = String("built in");
1220 }
1221//
1222 if (coeff.nelements()>0) {
1223 pPoly->setCoefficients(coeff);
1224 } else {
1225 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
1226 }
1227//
1228 cout << "Making polynomial correction with " << msg << " coefficients" << endl;
1229 const uInt nRow = in.nRow();
1230 Vector<Float> factor(nRow);
1231 for (uInt i=0; i<nRow; i++) {
1232 factor[i] = 1.0 / (*pPoly)(x[i]);
1233 }
1234 delete pPoly;
1235//
1236 scaleByVector (pTabOut, in, doAll, factor, True);
1237 } else {
1238
1239// Indicate which columns to read from ascii file
1240
1241 String col0("ELEVATION");
1242 String col1("FACTOR");
1243
1244// Read and correct
1245
1246 cout << "Making correction from ascii Table" << endl;
1247 scaleFromAsciiTable (pTabOut, in, fileName, col0, col1,
1248 methodStr, doAll, x, True);
1249 }
1250//
1251 return pTabOut;
1252}
1253
1254
1255
1256SDMemTable* SDMath::opacity(const SDMemTable& in, Float tau, Bool doAll) const
1257{
1258
1259// Get header and clone output table
1260
1261 SDHeader sh = in.getSDHeader();
1262 SDMemTable* pTabOut = new SDMemTable(in, True);
1263
1264// Get elevation data from SDMemTable and convert to degrees
1265
1266 const Table& tab = in.table();
1267 ROScalarColumn<Float> elev(tab, "ELEVATION");
1268 Vector<Float> zDist = elev.getColumn();
1269 zDist = Float(C::pi_2) - zDist;
1270
1271// Generate correction factor
1272
1273 const uInt nRow = in.nRow();
1274 Vector<Float> factor(nRow);
1275 Vector<Float> factor2(nRow);
1276 for (uInt i=0; i<nRow; i++) {
1277 factor[i] = exp(tau)/cos(zDist[i]);
1278 }
1279
1280// Correct
1281
1282 scaleByVector (pTabOut, in, doAll, factor, True);
1283//
1284 return pTabOut;
1285}
1286
1287
1288void SDMath::rotateXYPhase(SDMemTable& in, Float value, Bool doAll)
1289//
1290// phase in degrees
1291// assumes linear correlations
1292//
1293{
1294 if (in.nPol() != 4) {
1295 throw(AipsError("You must have 4 polarizations to run this function"));
1296 }
1297//
1298 SDHeader sh = in.getSDHeader();
1299 Instrument inst = SDAttr::convertInstrument (sh.antennaname, False);
1300 SDAttr sdAtt;
1301 if (sdAtt.feedPolType(inst) != LINEAR) {
1302 throw(AipsError("Only linear polarizations are supported"));
1303 }
1304//
1305 const Table& tabIn = in.table();
1306 ArrayColumn<Float> specCol(tabIn,"SPECTRA");
1307 IPosition start(asap::nAxes,0);
1308 IPosition end(asap::nAxes);
1309
1310// Set cursor slice. Assumes shape the same for all rows
1311
1312 setCursorSlice (start, end, doAll, in);
1313 IPosition start3(start);
1314 start3(asap::PolAxis) = 2; // Real(XY)
1315 IPosition end3(end);
1316 end3(asap::PolAxis) = 2;
1317//
1318 IPosition start4(start);
1319 start4(asap::PolAxis) = 3; // Imag (XY)
1320 IPosition end4(end);
1321 end4(asap::PolAxis) = 3;
1322//
1323 uInt nRow = in.nRow();
1324 Array<Float> data;
1325 for (uInt i=0; i<nRow;++i) {
1326 specCol.get(i,data);
1327 IPosition shape = data.shape();
1328
1329// Get polarization slice references
1330
1331 Array<Float> C3 = data(start3,end3);
1332 Array<Float> C4 = data(start4,end4);
1333
1334// Rotate
1335
1336 SDPolUtil::rotatePhase(C3, C4, value);
1337
1338// Put
1339
1340 specCol.put(i,data);
1341 }
1342}
1343
1344
1345
1346void SDMath::rotateLinPolPhase(SDMemTable& in, Float value, Bool doAll)
1347//
1348// phase in degrees
1349// assumes linear correlations
1350//
1351{
1352 if (in.nPol() != 4) {
1353 throw(AipsError("You must have 4 polarizations to run this function"));
1354 }
1355//
1356 SDHeader sh = in.getSDHeader();
1357 Instrument inst = SDAttr::convertInstrument (sh.antennaname, False);
1358 SDAttr sdAtt;
1359 if (sdAtt.feedPolType(inst) != LINEAR) {
1360 throw(AipsError("Only linear polarizations are supported"));
1361 }
1362//
1363 const Table& tabIn = in.table();
1364 ArrayColumn<Float> specCol(tabIn,"SPECTRA");
1365 ROArrayColumn<Float> stokesCol(tabIn,"STOKES");
1366 IPosition start(asap::nAxes,0);
1367 IPosition end(asap::nAxes);
1368
1369// Set cursor slice. Assumes shape the same for all rows
1370
1371 setCursorSlice (start, end, doAll, in);
1372//
1373 IPosition start1(start);
1374 start1(asap::PolAxis) = 0; // C1 (XX)
1375 IPosition end1(end);
1376 end1(asap::PolAxis) = 0;
1377//
1378 IPosition start2(start);
1379 start2(asap::PolAxis) = 1; // C2 (YY)
1380 IPosition end2(end);
1381 end2(asap::PolAxis) = 1;
1382//
1383 IPosition start3(start);
1384 start3(asap::PolAxis) = 2; // C3 ( Real(XY) )
1385 IPosition end3(end);
1386 end3(asap::PolAxis) = 2;
1387//
1388 IPosition startI(start);
1389 startI(asap::PolAxis) = 0; // I
1390 IPosition endI(end);
1391 endI(asap::PolAxis) = 0;
1392//
1393 IPosition startQ(start);
1394 startQ(asap::PolAxis) = 1; // Q
1395 IPosition endQ(end);
1396 endQ(asap::PolAxis) = 1;
1397//
1398 IPosition startU(start);
1399 startU(asap::PolAxis) = 2; // U
1400 IPosition endU(end);
1401 endU(asap::PolAxis) = 2;
1402
1403//
1404 uInt nRow = in.nRow();
1405 Array<Float> data, stokes;
1406 for (uInt i=0; i<nRow;++i) {
1407 specCol.get(i,data);
1408 stokesCol.get(i,stokes);
1409 IPosition shape = data.shape();
1410
1411// Get linear polarization slice references
1412
1413 Array<Float> C1 = data(start1,end1);
1414 Array<Float> C2 = data(start2,end2);
1415 Array<Float> C3 = data(start3,end3);
1416
1417// Get STokes slice references
1418
1419 Array<Float> I = stokes(startI,endI);
1420 Array<Float> Q = stokes(startQ,endQ);
1421 Array<Float> U = stokes(startU,endU);
1422
1423// Rotate
1424
1425 SDPolUtil::rotateLinPolPhase(C1, C2, C3, I, Q, U, value);
1426
1427// Put
1428
1429 specCol.put(i,data);
1430 }
1431}
1432
1433// 'private' functions
1434
1435void SDMath::convertBrightnessUnits (SDMemTable* pTabOut, const SDMemTable& in,
1436 Bool toKelvin, Float cFac, Bool doAll) const
1437{
1438
1439// Get header
1440
1441 SDHeader sh = in.getSDHeader();
1442 const uInt nChan = sh.nchan;
1443
1444// Get instrument
1445
1446 Bool throwIt = True;
1447 Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
1448
1449// Get Diameter (m)
1450
1451 SDAttr sdAtt;
1452
1453// Get epoch of first row
1454
1455 MEpoch dateObs = in.getEpoch(0);
1456
1457// Generate a Vector of correction factors. One per FreqID
1458
1459 SDFrequencyTable sdft = in.getSDFreqTable();
1460 Vector<uInt> freqIDs;
1461//
1462 Vector<Float> freqs(sdft.length());
1463 for (uInt i=0; i<sdft.length(); i++) {
1464 freqs(i) = (nChan/2 - sdft.referencePixel(i))*sdft.increment(i) + sdft.referenceValue(i);
1465 }
1466//
1467 Vector<Float> JyPerK = sdAtt.JyPerK(inst, dateObs, freqs);
1468 cout << "Jy/K = " << JyPerK << endl;
1469 Vector<Float> factors = cFac * JyPerK;
1470 if (toKelvin) factors = Float(1.0) / factors;
1471
1472// Get data slice bounds
1473
1474 IPosition start, end;
1475 setCursorSlice (start, end, doAll, in);
1476 const uInt ifAxis = in.getIF();
1477
1478// Iteration axes
1479
1480 IPosition axes(asap::nAxes-1,0);
1481 for (uInt i=0,j=0; i<asap::nAxes; i++) {
1482 if (i!=asap::IFAxis) {
1483 axes(j++) = i;
1484 }
1485 }
1486
1487// Loop over rows and apply correction factor
1488
1489 Float factor = 1.0;
1490 const uInt axis = asap::ChanAxis;
1491 for (uInt i=0; i < in.nRow(); ++i) {
1492
1493// Get data
1494
1495 MaskedArray<Float> dataIn = in.rowAsMaskedArray(i);
1496 Array<Float>& values = dataIn.getRWArray(); // Ref to dataIn
1497 Array<Float> values2 = values(start,end); // Ref to values to dataIn
1498
1499// Get SDCOntainer
1500
1501 SDContainer sc = in.getSDContainer(i);
1502
1503// Get FreqIDs
1504
1505 freqIDs = sc.getFreqMap();
1506
1507// Now the conversion factor depends only upon frequency
1508// So we need to iterate through by IF only giving
1509// us BEAM/POL/CHAN cubes
1510
1511 ArrayIterator<Float> itIn(values2, axes);
1512 uInt ax = 0;
1513 while (!itIn.pastEnd()) {
1514 itIn.array() *= factors(freqIDs(ax)); // Writes back to dataIn
1515 itIn.next();
1516 }
1517
1518// Write out
1519
1520 putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
1521//
1522 pTabOut->putSDContainer(sc);
1523 }
1524}
1525
1526
1527
1528SDMemTable* SDMath::frequencyAlign (const SDMemTable& in,
1529 MFrequency::Types freqSystem,
1530 const String& refTime,
1531 const String& methodStr,
1532 Bool perFreqID) const
1533{
1534// Get Header
1535
1536 SDHeader sh = in.getSDHeader();
1537 const uInt nChan = sh.nchan;
1538 const uInt nRows = in.nRow();
1539 const uInt nIF = sh.nif;
1540
1541// Get Table reference
1542
1543 const Table& tabIn = in.table();
1544
1545// Get Columns from Table
1546
1547 ROScalarColumn<Double> mjdCol(tabIn, "TIME");
1548 ROScalarColumn<String> srcCol(tabIn, "SRCNAME");
1549 ROArrayColumn<uInt> fqIDCol(tabIn, "FREQID");
1550 Vector<Double> times = mjdCol.getColumn();
1551
1552// Generate DataDesc table
1553
1554 Matrix<uInt> ddIdx;
1555 SDDataDesc dDesc;
1556 generateDataDescTable (ddIdx, dDesc, nIF, in, tabIn, srcCol, fqIDCol, perFreqID);
1557
1558// Get reference Epoch to time of first row or given String
1559
1560 Unit DAY(String("d"));
1561 MEpoch::Ref epochRef(in.getTimeReference());
1562 MEpoch refEpoch;
1563 if (refTime.length()>0) {
1564 refEpoch = epochFromString(refTime, in.getTimeReference());
1565 } else {
1566 refEpoch = in.getEpoch(0);
1567 }
1568 cout << "Aligning at reference Epoch " << formatEpoch(refEpoch)
1569 << " in frame " << MFrequency::showType(freqSystem) << endl;
1570
1571// Get Reference Position
1572
1573 MPosition refPos = in.getAntennaPosition();
1574
1575// Create FrequencyAligner Block. One FA for each possible
1576// source/freqID (perFreqID=True) or source/IF (perFreqID=False) combination
1577
1578 PtrBlock<FrequencyAligner<Float>* > a(dDesc.length());
1579 generateFrequencyAligners (a, dDesc, in, nChan, freqSystem, refPos,
1580 refEpoch, perFreqID);
1581
1582// Generate and fill output Frequency Table. WHen perFreqID=True, there is one output FreqID
1583// for each entry in the SDDataDesc table. However, in perFreqID=False mode, there may be
1584// some degeneracy, so we need a little translation map
1585
1586 SDFrequencyTable freqTabOut = in.getSDFreqTable();
1587 freqTabOut.setLength(0);
1588 Vector<String> units(1);
1589 units = String("Hz");
1590 Bool linear=True;
1591//
1592 Vector<uInt> ddFQTrans(dDesc.length(),0);
1593 for (uInt i=0; i<dDesc.length(); i++) {
1594
1595// Get Aligned SC in Hz
1596
1597 SpectralCoordinate sC = a[i]->alignedSpectralCoordinate(linear);
1598 sC.setWorldAxisUnits(units);
1599
1600// Add FreqID
1601
1602 uInt idx = freqTabOut.addFrequency(sC.referencePixel()[0],
1603 sC.referenceValue()[0],
1604 sC.increment()[0]);
1605 ddFQTrans(i) = idx; // output FreqID = ddFQTrans(ddIdx)
1606 }
1607
1608// Interpolation method
1609
1610 InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1611 convertInterpString(interp, methodStr);
1612
1613// New output Table
1614
1615 cout << "Create output table" << endl;
1616 SDMemTable* pTabOut = new SDMemTable(in,True);
1617 pTabOut->putSDFreqTable(freqTabOut);
1618
1619// Loop over rows in Table
1620
1621 Bool extrapolate=False;
1622 const IPosition polChanAxes(2, asap::PolAxis, asap::ChanAxis);
1623 Bool useCachedAbcissa = False;
1624 Bool first = True;
1625 Bool ok;
1626 Vector<Float> yOut;
1627 Vector<Bool> maskOut;
1628 Vector<uInt> freqID(nIF);
1629 uInt ifIdx, faIdx;
1630 Vector<Double> xIn;
1631//
1632 for (uInt iRow=0; iRow<nRows; ++iRow) {
1633 if (iRow%10==0) {
1634 cout << "Processing row " << iRow << endl;
1635 }
1636
1637// Get EPoch
1638
1639 Quantum<Double> tQ2(times[iRow],DAY);
1640 MVEpoch mv2(tQ2);
1641 MEpoch epoch(mv2, epochRef);
1642
1643// Get copy of data
1644
1645 const MaskedArray<Float>& mArrIn(in.rowAsMaskedArray(iRow));
1646 Array<Float> values = mArrIn.getArray();
1647 Array<Bool> mask = mArrIn.getMask();
1648
1649// For each row, the Frequency abcissa will be the same regardless
1650// of polarization. For all other axes (IF and BEAM) the abcissa
1651// will change. So we iterate through the data by pol-chan planes
1652// to mimimize the work. Probably won't work for multiple beams
1653// at this point.
1654
1655 ArrayIterator<Float> itValuesPlane(values, polChanAxes);
1656 ArrayIterator<Bool> itMaskPlane(mask, polChanAxes);
1657 while (!itValuesPlane.pastEnd()) {
1658
1659// Find the IF index and then the FA PtrBlock index
1660
1661 const IPosition& pos = itValuesPlane.pos();
1662 ifIdx = pos(asap::IFAxis);
1663 faIdx = ddIdx(iRow,ifIdx);
1664
1665// Generate abcissa for perIF. Could cache this in a Matrix
1666// on a per scan basis. Pretty expensive doing it for every row.
1667
1668 if (!perFreqID) {
1669 xIn.resize(nChan);
1670 uInt fqID = dDesc.secID(ddIdx(iRow,ifIdx));
1671 SpectralCoordinate sC = in.getSpectralCoordinate(fqID);
1672 Double w;
1673 for (uInt i=0; i<nChan; i++) {
1674 sC.toWorld(w,Double(i));
1675 xIn[i] = w;
1676 }
1677 }
1678//
1679 VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1);
1680 VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
1681
1682// Iterate through the plane by vector and align
1683
1684 first = True;
1685 useCachedAbcissa=False;
1686 while (!itValuesVec.pastEnd()) {
1687 if (perFreqID) {
1688 ok = a[faIdx]->align (yOut, maskOut, itValuesVec.vector(),
1689 itMaskVec.vector(), epoch, useCachedAbcissa,
1690 interp, extrapolate);
1691 } else {
1692 ok = a[faIdx]->align (yOut, maskOut, xIn, itValuesVec.vector(),
1693 itMaskVec.vector(), epoch, useCachedAbcissa,
1694 interp, extrapolate);
1695 }
1696//
1697 itValuesVec.vector() = yOut;
1698 itMaskVec.vector() = maskOut;
1699//
1700 itValuesVec.next();
1701 itMaskVec.next();
1702//
1703 if (first) {
1704 useCachedAbcissa = True;
1705 first = False;
1706 }
1707 }
1708//
1709 itValuesPlane.next();
1710 itMaskPlane.next();
1711 }
1712
1713// Create SDContainer and put back
1714
1715 SDContainer sc = in.getSDContainer(iRow);
1716 putDataInSDC(sc, values, mask);
1717
1718// Set output FreqIDs
1719
1720 for (uInt i=0; i<nIF; i++) {
1721 uInt idx = ddIdx(iRow,i); // Index into SDDataDesc table
1722 freqID(i) = ddFQTrans(idx); // FreqID in output FQ table
1723 }
1724 sc.putFreqMap(freqID);
1725//
1726 pTabOut->putSDContainer(sc);
1727 }
1728
1729// Now we must set the base and extra frames to the
1730// input frame
1731
1732 std::vector<string> info = pTabOut->getCoordInfo();
1733 info[1] = MFrequency::showType(freqSystem); // Conversion frame
1734 info[3] = info[1]; // Base frame
1735 pTabOut->setCoordInfo(info);
1736
1737// Clean up PointerBlock
1738
1739 for (uInt i=0; i<a.nelements(); i++) delete a[i];
1740//
1741 return pTabOut;
1742}
1743
1744
1745void SDMath::fillSDC(SDContainer& sc,
1746 const Array<Bool>& mask,
1747 const Array<Float>& data,
1748 const Array<Float>& tSys,
1749 Int scanID, Double timeStamp,
1750 Double interval, const String& sourceName,
1751 const Vector<uInt>& freqID) const
1752{
1753// Data and mask
1754
1755 putDataInSDC(sc, data, mask);
1756
1757// TSys
1758
1759 sc.putTsys(tSys);
1760
1761// Time things
1762
1763 sc.timestamp = timeStamp;
1764 sc.interval = interval;
1765 sc.scanid = scanID;
1766//
1767 sc.sourcename = sourceName;
1768 sc.putFreqMap(freqID);
1769}
1770
1771void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum,
1772 MaskedArray<Float>& sum, Array<Float>& sumSq,
1773 Array<Float>& nPts, Array<Float>& tSysSum,
1774 Array<Float>& tSysSqSum,
1775 const Array<Float>& tSys, const Array<Float>& nInc,
1776 const Vector<Bool>& mask, Double time, Double interval,
1777 const Block<CountedPtr<SDMemTable> >& in,
1778 uInt iTab, uInt iRow, uInt axis,
1779 uInt nAxesSub, Bool useMask,
1780 WeightType wtType) const
1781{
1782
1783// Get data
1784
1785 MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
1786 Array<Float>& valuesIn = dataIn.getRWArray(); // writable reference
1787 const Array<Bool>& maskIn = dataIn.getMask(); // RO reference
1788//
1789 if (wtType==NONE) {
1790 const MaskedArray<Float> n(nInc,dataIn.getMask());
1791 nPts += n; // Only accumulates where mask==T
1792 } else if (wtType==TINT) {
1793
1794// We are weighting the data by integration time.
1795
1796 valuesIn *= Float(interval);
1797
1798 } else if (wtType==VAR) {
1799
1800// We are going to average the data, weighted by the noise for each pol, beam and IF.
1801// So therefore we need to iterate through by spectrum (axis 3)
1802
1803 VectorIterator<Float> itData(valuesIn, axis);
1804 ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
1805 Float fac = 1.0;
1806 IPosition pos(nAxesSub,0);
1807//
1808 while (!itData.pastEnd()) {
1809
1810// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
1811
1812 if (useMask) {
1813 MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
1814 fac = 1.0/variance(tmp);
1815 } else {
1816 MaskedArray<Float> tmp(itData.vector(),itMask.vector());
1817 fac = 1.0/variance(tmp);
1818 }
1819
1820// Scale data
1821
1822 itData.vector() *= fac; // Writes back into 'dataIn'
1823//
1824// Accumulate variance per if/pol/beam averaged over spectrum
1825// This method to get pos2 from itData.pos() is only valid
1826// because the spectral axis is the last one (so we can just
1827// copy the first nAXesSub positions out)
1828
1829 pos = itData.pos().getFirst(nAxesSub);
1830 sumSq(pos) += fac;
1831//
1832 itData.next();
1833 itMask.next();
1834 }
1835 } else if (wtType==TSYS) {
1836
1837// We are going to average the data, weighted by 1/Tsys**2 for each pol, beam and IF.
1838// So therefore we need to iterate through by spectrum (axis 3). Although
1839// Tsys is stored as a vector of length nChan, the values are replicated.
1840// We will take a short cut and just use the value from the first channel
1841// for now.
1842//
1843 VectorIterator<Float> itData(valuesIn, axis);
1844 ReadOnlyVectorIterator<Float> itTSys(tSys, axis);
1845 Float fac = 1.0;
1846 IPosition pos(nAxesSub,0);
1847//
1848 while (!itData.pastEnd()) {
1849 Float t = itTSys.vector()[0];
1850 fac = 1.0/t/t;
1851
1852// Scale data
1853
1854 itData.vector() *= fac; // Writes back into 'dataIn'
1855//
1856// Accumulate Tsys per if/pol/beam averaged over spectrum
1857// This method to get pos2 from itData.pos() is only valid
1858// because the spectral axis is the last one (so we can just
1859// copy the first nAXesSub positions out)
1860
1861 pos = itData.pos().getFirst(nAxesSub);
1862 tSysSqSum(pos) += fac;
1863//
1864 itData.next();
1865 itTSys.next();
1866 }
1867 }
1868
1869// Accumulate sum of (possibly scaled) data
1870
1871 sum += dataIn;
1872
1873// Accumulate Tsys, time, and interval
1874
1875 tSysSum += tSys;
1876 timeSum += time;
1877 intSum += interval;
1878 nAccum += 1;
1879}
1880
1881
1882void SDMath::normalize(MaskedArray<Float>& sum,
1883 const Array<Float>& sumSq,
1884 const Array<Float>& tSysSqSum,
1885 const Array<Float>& nPts,
1886 Double intSum,
1887 WeightType wtType, Int axis,
1888 Int nAxesSub) const
1889{
1890 IPosition pos2(nAxesSub,0);
1891//
1892 if (wtType==NONE) {
1893
1894// We just average by the number of points accumulated.
1895// We need to make a MA out of nPts so that no divide by
1896// zeros occur
1897
1898 MaskedArray<Float> t(nPts, (nPts>Float(0.0)));
1899 sum /= t;
1900 } else if (wtType==TINT) {
1901
1902// Average by sum of Tint
1903
1904 sum /= Float(intSum);
1905 } else if (wtType==VAR) {
1906
1907// Normalize each spectrum by sum(1/var) where the variance
1908// is worked out for each spectrum
1909
1910 Array<Float>& data = sum.getRWArray();
1911 VectorIterator<Float> itData(data, axis);
1912 while (!itData.pastEnd()) {
1913 pos2 = itData.pos().getFirst(nAxesSub);
1914 itData.vector() /= sumSq(pos2);
1915 itData.next();
1916 }
1917 } else if (wtType==TSYS) {
1918
1919// Normalize each spectrum by sum(1/Tsys**2) where the pseudo
1920// replication over channel for Tsys has been dropped.
1921
1922 Array<Float>& data = sum.getRWArray();
1923 VectorIterator<Float> itData(data, axis);
1924 while (!itData.pastEnd()) {
1925 pos2 = itData.pos().getFirst(nAxesSub);
1926 itData.vector() /= tSysSqSum(pos2);
1927 itData.next();
1928 }
1929 }
1930}
1931
1932
1933
1934
1935void SDMath::setCursorSlice (IPosition& start, IPosition& end, Bool doAll, const SDMemTable& in) const
1936{
1937 const uInt nDim = asap::nAxes;
1938 DebugAssert(nDim==4,AipsError);
1939//
1940 start.resize(nDim);
1941 end.resize(nDim);
1942 if (doAll) {
1943 start = 0;
1944 end(0) = in.nBeam()-1;
1945 end(1) = in.nIF()-1;
1946 end(2) = in.nPol()-1;
1947 end(3) = in.nChan()-1;
1948 } else {
1949 start(0) = in.getBeam();
1950 end(0) = start(0);
1951//
1952 start(1) = in.getIF();
1953 end(1) = start(1);
1954//
1955 start(2) = in.getPol();
1956 end(2) = start(2);
1957//
1958 start(3) = 0;
1959 end(3) = in.nChan()-1;
1960 }
1961}
1962
1963
1964void SDMath::convertWeightString(WeightType& wtType, const String& weightStr,
1965 Bool listType) const
1966{
1967 String tStr(weightStr);
1968 tStr.upcase();
1969 String msg;
1970 if (tStr.contains(String("NONE"))) {
1971 wtType = NONE;
1972 msg = String("Weighting type selected : None");
1973 } else if (tStr.contains(String("VAR"))) {
1974 wtType = VAR;
1975 msg = String("Weighting type selected : Variance");
1976 } else if (tStr.contains(String("TINT"))) {
1977 wtType = TINT;
1978 msg = String("Weighting type selected : Tint");
1979 } else if (tStr.contains(String("TSYS"))) {
1980 wtType = TSYS;
1981 msg = String("Weighting type selected : Tsys");
1982 } else {
1983 msg = String("Weighting type selected : None");
1984 throw(AipsError("Unrecognized weighting type"));
1985 }
1986//
1987 if (listType) cout << msg << endl;
1988}
1989
1990
1991void SDMath::convertInterpString(casa::InterpolateArray1D<Double,Float>::InterpolationMethod& type,
1992 const casa::String& interp) const
1993{
1994 String tStr(interp);
1995 tStr.upcase();
1996 if (tStr.contains(String("NEAR"))) {
1997 type = InterpolateArray1D<Double,Float>::nearestNeighbour;
1998 } else if (tStr.contains(String("LIN"))) {
1999 type = InterpolateArray1D<Double,Float>::linear;
2000 } else if (tStr.contains(String("CUB"))) {
2001 type = InterpolateArray1D<Double,Float>::cubic;
2002 } else if (tStr.contains(String("SPL"))) {
2003 type = InterpolateArray1D<Double,Float>::spline;
2004 } else {
2005 throw(AipsError("Unrecognized interpolation type"));
2006 }
2007}
2008
2009void SDMath::putDataInSDC(SDContainer& sc, const Array<Float>& data,
2010 const Array<Bool>& mask) const
2011{
2012 sc.putSpectrum(data);
2013//
2014 Array<uChar> outflags(data.shape());
2015 convertArray(outflags,!mask);
2016 sc.putFlags(outflags);
2017}
2018
2019Table SDMath::readAsciiFile (const String& fileName) const
2020{
2021 String formatString;
2022 Table tbl = readAsciiTable (formatString, Table::Memory, fileName, "", "", False);
2023 return tbl;
2024}
2025
2026
2027
2028void SDMath::scaleFromAsciiTable(SDMemTable* pTabOut,
2029 const SDMemTable& in, const String& fileName,
2030 const String& col0, const String& col1,
2031 const String& methodStr, Bool doAll,
2032 const Vector<Float>& xOut, Bool doTSys) const
2033{
2034
2035// Read gain-elevation ascii file data into a Table.
2036
2037 Table geTable = readAsciiFile (fileName);
2038//
2039 scaleFromTable (pTabOut, in, geTable, col0, col1, methodStr, doAll, xOut, doTSys);
2040}
2041
2042void SDMath::scaleFromTable(SDMemTable* pTabOut, const SDMemTable& in,
2043 const Table& tTable, const String& col0,
2044 const String& col1,
2045 const String& methodStr, Bool doAll,
2046 const Vector<Float>& xOut, Bool doTsys) const
2047{
2048
2049// Get data from Table
2050
2051 ROScalarColumn<Float> geElCol(tTable, col0);
2052 ROScalarColumn<Float> geFacCol(tTable, col1);
2053 Vector<Float> xIn = geElCol.getColumn();
2054 Vector<Float> yIn = geFacCol.getColumn();
2055 Vector<Bool> maskIn(xIn.nelements(),True);
2056
2057// Interpolate (and extrapolate) with desired method
2058
2059 InterpolateArray1D<Double,Float>::InterpolationMethod method;
2060 convertInterpString(method, methodStr);
2061 Int intMethod(method);
2062//
2063 Vector<Float> yOut;
2064 Vector<Bool> maskOut;
2065 InterpolateArray1D<Float,Float>::interpolate(yOut, maskOut, xOut,
2066 xIn, yIn, maskIn, intMethod,
2067 True, True);
2068// Apply
2069
2070 scaleByVector(pTabOut, in, doAll, Float(1.0)/yOut, doTsys);
2071}
2072
2073
2074void SDMath::scaleByVector(SDMemTable* pTabOut, const SDMemTable& in,
2075 Bool doAll, const Vector<Float>& factor,
2076 Bool doTSys) const
2077{
2078
2079// Set up data slice
2080
2081 IPosition start, end;
2082 setCursorSlice (start, end, doAll, in);
2083
2084// Get Tsys column
2085
2086 const Table& tIn = in.table();
2087 ArrayColumn<Float> tSysCol(tIn, "TSYS");
2088 Array<Float> tSys;
2089
2090// Loop over rows and apply correction factor
2091
2092 const uInt axis = asap::ChanAxis;
2093 for (uInt i=0; i < in.nRow(); ++i) {
2094
2095// Get data
2096
2097 MaskedArray<Float> dataIn(in.rowAsMaskedArray(i));
2098 MaskedArray<Float> dataIn2 = dataIn(start,end); // reference to dataIn
2099//
2100 if (doTSys) {
2101 tSysCol.get(i, tSys);
2102 Array<Float> tSys2 = tSys(start,end) * factor[i];
2103 tSysCol.put(i, tSys);
2104 }
2105
2106// Apply factor
2107
2108 dataIn2 *= factor[i];
2109
2110// Write out
2111
2112 SDContainer sc = in.getSDContainer(i);
2113 putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
2114//
2115 pTabOut->putSDContainer(sc);
2116 }
2117}
2118
2119
2120
2121
2122void SDMath::generateDataDescTable (Matrix<uInt>& ddIdx,
2123 SDDataDesc& dDesc,
2124 uInt nIF,
2125 const SDMemTable& in,
2126 const Table& tabIn,
2127 const ROScalarColumn<String>& srcCol,
2128 const ROArrayColumn<uInt>& fqIDCol,
2129 Bool perFreqID) const
2130{
2131 const uInt nRows = tabIn.nrow();
2132 ddIdx.resize(nRows,nIF);
2133//
2134 String srcName;
2135 Vector<uInt> freqIDs;
2136 for (uInt iRow=0; iRow<nRows; iRow++) {
2137 srcCol.get(iRow, srcName);
2138 fqIDCol.get(iRow, freqIDs);
2139 const MDirection& dir = in.getDirection(iRow);
2140//
2141 if (perFreqID) {
2142
2143// One entry per source/freqID pair
2144
2145 for (uInt iIF=0; iIF<nIF; iIF++) {
2146 ddIdx(iRow,iIF) = dDesc.addEntry(srcName, freqIDs[iIF], dir, 0);
2147 }
2148 } else {
2149
2150// One entry per source/IF pair. Hang onto the FreqID as well
2151
2152 for (uInt iIF=0; iIF<nIF; iIF++) {
2153 ddIdx(iRow,iIF) = dDesc.addEntry(srcName, iIF, dir, freqIDs[iIF]);
2154 }
2155 }
2156 }
2157}
2158
2159
2160
2161
2162
2163MEpoch SDMath::epochFromString (const String& str, MEpoch::Types timeRef) const
2164{
2165 Quantum<Double> qt;
2166 if (MVTime::read(qt,str)) {
2167 MVEpoch mv(qt);
2168 MEpoch me(mv, timeRef);
2169 return me;
2170 } else {
2171 throw(AipsError("Invalid format for Epoch string"));
2172 }
2173}
2174
2175
2176String SDMath::formatEpoch(const MEpoch& epoch) const
2177{
2178 MVTime mvt(epoch.getValue());
2179 return mvt.string(MVTime::YMD) + String(" (") + epoch.getRefString() + String(")");
2180}
2181
2182
2183
2184void SDMath::generateFrequencyAligners (PtrBlock<FrequencyAligner<Float>* >& a,
2185 const SDDataDesc& dDesc,
2186 const SDMemTable& in, uInt nChan,
2187 MFrequency::Types system,
2188 const MPosition& refPos,
2189 const MEpoch& refEpoch,
2190 Bool perFreqID) const
2191{
2192 for (uInt i=0; i<dDesc.length(); i++) {
2193 uInt ID = dDesc.ID(i);
2194 uInt secID = dDesc.secID(i);
2195 const MDirection& refDir = dDesc.secDir(i);
2196//
2197 if (perFreqID) {
2198
2199// One aligner per source/FreqID pair.
2200
2201 SpectralCoordinate sC = in.getSpectralCoordinate(ID);
2202 a[i] = new FrequencyAligner<Float>(sC, nChan, refEpoch, refDir, refPos, system);
2203 } else {
2204
2205// One aligner per source/IF pair. But we still need the FreqID to
2206// get the right SC. Hence the messing about with the secondary ID
2207
2208 SpectralCoordinate sC = in.getSpectralCoordinate(secID);
2209 a[i] = new FrequencyAligner<Float>(sC, nChan, refEpoch, refDir, refPos, system);
2210 }
2211 }
2212}
2213
2214Vector<uInt> SDMath::getRowRange (const SDMemTable& in) const
2215{
2216 Vector<uInt> range(2);
2217 range[0] = 0;
2218 range[1] = in.nRow()-1;
2219 return range;
2220}
2221
2222
2223Bool SDMath::rowInRange (uInt i, const Vector<uInt>& range) const
2224{
2225 return (i>=range[0] && i<=range[1]);
2226}
Note: See TracBrowser for help on using the repository browser.