source: tags/Release-2/src/SDMath.cc@ 1439

Last change on this file since 1439 was 653, checked in by mar637, 19 years ago

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