source: branches/Release-1-fixes/src/SDMath.cc@ 2468

Last change on this file since 2468 was 536, checked in by kil064, 20 years ago

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